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ABSTRACT 


The  Diffusive  Transport  Module  of  the  new  DELFIC  fall¬ 
out  prediction  system  has  undergone  additional  development 
since  publication  of  its  description  in  DASA  2669.  This 
supplement  to  DASA  2669  describes  these  developments  and 
presents  ammendments  and  corrections  to  the  code  ar.d  its 
documentation.  Complete  FORTRAN  statement  listings  of  sub¬ 
routines  that  have  been  substantially  changed  are  included. 
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1.  INTRODUCTION 

Since  publication  of  DASA  2669*,  development  and  appli¬ 
cation  of  the  DELFIC  Diffusive  Transport  Module  (DTM)  has  continued. 

A  few  revisions  to  the  nodel  have  been  made,  several  of  them  of 
major  importance.  In  addition,  some  hidden  "bugs"  have  been  uncovered 
and  corrected.  In  this  supplement  to  DASA  2669  we  describe  the 
important  model  revisions,  and  amend  the  documentation.  We  alco 
correct  errors  in  the  documentation,  and  provide  FORTRAN  statement 
listings  of  subroutines  that  have  been  changed  substantially. 


*  H.  G.  Norment  ana  E.  J.  Tichovolsky,  "A  New  Fallout  Transport  Code 
for  the  DELFIC  System:  The  Diffusive  Transport  Module,"  ARCON 
Corporation  Report  R71-1W,  DASA  2669  (1  March  1973),  AD  727  613. 


1 


»«“*«'<»  itfjy  h 


2.  MODEL  REVISIONS 


2.1  Initial  Parcel  Description 

Fallout  parcels  are  taken  to  be  distributed  in  the  hori¬ 
zontal  about  their  centers  of  mass  by  a  Gaussian  density  function. 

The  initial  Gaussian  standard  deviation  was  set  equal  to  the  parcel 
radius  that  was  received  via  tape  IPASIN  from  the  Cloud  Rise-Transport 
Interface  Module.  This  has  been  changed  so  that  the  initial  standard 
deviation  is  one-half  the  input  parcel  radius.  (Compare  subroutine 
SPRVS  card  81  of  DASA  2669  with  SPRVS  card  94  of  this  supplement.) 


2.2  Simple  Advection-Plus-Settling 

In  some  cases,  it  is  desirable  to  transport  fallout  in  a 
simple  advection-plus-settling  mode;  that  is,  without  accounting  for 
diffusion  in  the  vertical.  In  this  mode,  integration  of  Eq.  (16)  is 
bypassed,  and  the  parcel  trajectory  is  computed  via  Eq.  (32).  As 
actually  employed  in  the  original  DTM,  Eq.  (32)  was  modified  to  the 
form 

z 

g 

i  +  <f>  l  <w>  ^  ^(x’y>2't)Az  ’ 

z . 

i 

where  the  average  settling  speed,  <f>,  was  taken  to  be 


<f>  = 


f(z.)  +  f (z  ) 


and  <w>  was  an  average  vertical  air  velocity.  For  cases  where  z.-z 
is  large,  there  can  be  significant  differences  in  the  particle 
settling  speeds  in  the  upper  layers  compared  with  those  in  the  lower 
layers.  When,  in  addition,  there  is  large  wind  shear  in  the  vertical. 


Preceding  page  blank 
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it  becomes  necessary  to  use  settling  speeds  that  are  computed  indi¬ 
vidually  for  p:Ji  -ind  layer  instead  of  an  average  settling  speed. 

The  code  has  been  changed  so  that  this  is  done. 

Major  changes  have  been  required  in  several  programs.  The 
most  extensive  changes  are  in  subroutines  SPRVS  and  TRANP.  A  new 
array,  CAVS(KBHF) ,  is  created  to  store  a  table  of  particle  settling 
speeds;  it  contains  an  entry  for  each  wind  layer. 

Before  a  parcel  is  transported  vir  this  mode,  a  test  is 
made  to  determine  if  the  parcel  can  impact  in  the  time  allowed  for 
transport.  This  test  is  simplified  by  using  the  knowledge  that  all 
parcels  comprised  of  particles  of  a  particular  size  class  are  processed 
sequentially  in  one  group.  Thus,  as  each  new  particle  size  class  is 
encountered,  the  a^ticude  above  which  these  particles  cannot  impact  is 
computed.  Then,  any  parcel  in  the  group  whose  base  is  above  cnis  alti¬ 
tude  is  bypassed.  To  perform  the  altitude  limit  calculation,  an 
average  vertical  wind  velocity  is  needed  for  each  wind  layer.  To 
accommodate  this,  the  array  WAVG(LTIMF)  was  changed  to  the  array 
WAVG  (  KBHK ,  I»1  IMF )  . 


3.  DOCUMENTATION  REVISIONS  ; 

Revisions  and  corrections  are  intermixed  and  listed  in 
order  of  their  encounter  in  DASA  2669. 


Page  27,  Eq.  (23): 


£  =  4c 

C1  3 


t. 

i 


1/3 


(23) 


Page  27,  line  16: 

In  our  application,  ct  is  taken  to  be  one-half  of  the 

i 

radius  of  a  cloud  wafer  as  recorded 


Page  44,  line  3: 

where  the  summation  is  over  the  N  data  with  the  largest  f^, 
and  the  weighting  factors, 


m  \ 


Page  44,  lines  6  and  1: 

The  parameters  a,  (3,  and  N  are  specified  by  the  user  and 

the  x.  and  y.  are  relative  to  the  n-th  lattice  cell  center, 
l  i 

The  calculations  of  the  f  are  performed  so  that  whenever  a 
factor  in  Eq.  (39)  is  found  to 

Page  44,  lines  9  and  10: 

than  the  total  number  of  observations,  M,  only  the  N  obser¬ 
vations  with  the  largest  f  relative  to  the  n-th  lattice  cell 
center  are  considered  in  the  calculations. 

Page  45,  lines  8-12: 

list,  two  tables  of  settling  rates  for  this  particle  are 
computed.  One  table  contains  an  entry  for  each  altitude 
increment  used  in  the  numerical  integration  of  Eq.  (16). 
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The  other  table  contains  an  entry  for  each  of  the  larger 
altitude  increments  that  are  used  to  resolve  the  wind  and 
turbulence  fields.  These  same  tables  are  used  through 
transport  for  this  parcel  and  the  remainder  of  the  parcels 
in  the  first  group.  When  the  first  parcel  of  the  second 
group  is  selected,  new  tables  are  computed,  and  so  on. 


Page  49,  Eq.  (45): 


Gn+1  -  G" 

J  J  _  1 


At 


2(Ae 


„\ 2  !(Kj+i +  Wi 


(K.  ,  +  2K.  +  K.  ,)Gn 
J+l  J  3-1  J 


(K.  +  K.  .)Gn  .I+7-  [(f.^.  -  w.  .)Gn  -  (f.  -  w . ) G^l  (45) 
3  J-l  j-lJ  Az  l_v  j+l  j+l  j+l  j  J  Jj 


Page  53,  following  Eq.  (65): 


where  in  the  code  is  taken  to  be  equal  to  K^. 


Page  66,  paragraphs  1  and  2: 

Beginning  at  its  input  location  and  time,  the  parcel 
base  or  top  is  transported  via  local  winds  in  the  cell  of  its 
residence.  At  the  same  time  it  settles  at  a  speed  that  is 
computed  for  the  altitude  at  the  center  of  the  wind  cell  of 
its  residence.  When  it  passes  through  a  wind  cell  boundary 
or  a  time  boundary,  the  wind  and  settling  speed  are  changed 
to  those  of  the  new  cell  or  update.  This  continues  until 
ground  impaction  occurs  or  until  an  extreme  wind  fieid  or 
time  boundary  is  encountered.  The  calculation  requires  one 
step  for  each  cell  through  which  the  base  or  top  passes. 

When  an  extreme  boundary  is  encountered  by  a  top  or  base, 
the  location  and  time  of  the  encounter  is  recorded.  These 
values  are  used  in  the  definition  of  the  deposit  increment  ?s 
described  on  page  70.  The  altitude  of  a  deposit  increment  is 
always  recorded  as  the  arithmetic  average  of  the  impact 


altitudes  of  its  top  and  base.  Thus,  the  recorded  alti¬ 
tude  of  a  deposit  increment  that  has  reached  an  extreme 
boundary  can  be  veil  above  the  deposition  plane. 

Page  71,  Figure  9: 

The  quantity  labeled  cr(J_l_  that  lies  to  the  left  of  the 
deposit  increment  ellipse  should  be  replaced  by  °r(_[_)c)* 

Page  88,  lines  10  and  11: 

DFZ.  An  area-weighted  average  vertical  wind,  WAVG(KBH,  LTIM) , 
is  derived  from  array  WFZ  for  each  altitude  layer  and  update. 
Likewise,  a  volume -weigh ted 

Page  88,  equation  for  DKAV: 


DKAV  = 


KBHX-1 

1 


KBH=1 


DFZ (KBH.NDATA. LTIM) *(ZBH(KBml)  -  ZBH(KBH)) 
ZBH(KBHX)  -  ZBH(l) 


Page  88,  lines  27  and  28: 

out  on  ISOUT.  In  a  parallel  operation  the  quantities 
WAVG(KBH,LTIM) ,  which  are  area-weighted  vertical  wind 
velocities  for  each  wind  layer  and  update,  are  computed 
and  printed  out  on 


Page  96,  line  25: 

in  the  horizontal  RWFR(J)/2.0; 


page  9V,  line  12: 

ZPAR(J) ,  ?SAM(J),  RWFR(J) /2 .0,  DWFR(J) ,  ZLWF(J),  and 
VWFR(J)  are 


Page  98,  line  2: 

for  further  details.)  Also  computed  are  the  settling  speeds. 
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CAVS(KBH),  for  each  altitude,  and  the  altitude,  ZLIM,  above 
which  deposition  is  impossible,  via  gravity  settling  in  the 
specified  wind  field,  for  particles  in  the  particular  size 
class  being  considered. 

Page  100,  lasc  line: 

CAY  -  CAV  -  K'AVGK 


Page  101,  first  line: 

WAVGK  and  DAVG^l)  are  the  average  vertical  components  of  wind 
Page  101,  line  15: 

ground  by  advection  at  fall  rates  CAVS(KBHF)  via  a  call  to 
subroutine  ADVEC. 

i 

Page  101,  line  23: 

limit:  i.e.,  whenever  ZLCW  >  ZLIM.  The  comment 

Page  110,  lines  16-19:  ’ 

also  obtained  via  the  COMMON  area  QPARM.  The  particle 
settling  speeds,  CAVS(KBHF),  and  area-weighted  vertical 
air  velocities,  WAYG(KBHT,LTIMF) ,  for  each- wind  layer,  are 
obtained 

Page  110,  line  27: 

parcel  base  is  advected,  while  settling  at  speeds  CAVS(KBH)  - 
WAVG (KBH, LTIM) ,  from  position  (XP,YP,ZP  = 


Page  lip,  line  29: 

TOL.  The  standard  deviations  of  the 


8 


Page  110,  line  33: 

parcel  top  if  advected,  while  settling  at  speeds  CAVSlKBH;  - 
WAVG(KBH,LTI?0 ,  from  position  (XP,Yr, 

Page  111,  line  2: 

at  time  TOU.  The  standard 

Page  111,  second  paragraph: 

After  both  base  and  top  of  a  parcel  have  been  transported, 
the  arithmetic  average  of  the  two  impact  times  and  impact 
altitudes  are  recorded  for  the  deposit  increment  (as  show-* 
on  the  next  page).  In  the  event  that  the  base  or  top  encounters 
an  extreme  wind  field  or  time  boundary  during  transport,  sub- 
soutine  TRANP  returns  the  coordinates  of  the  encounter  point. 
Therefore,  the  altitude  recorded  for  a  deposit  increment  can 
be  well  above  the  deposition  plane. 

Page  111,  third  paragraph,  lines  4  and  5: 

is  considered  superfluous  whenever  ZP  differs  from  ZL-cP  by 
less  than  0.1.  Instead,  XOL,  YOL,  ZOL,  TGL,  SIGXT-,  SIGYL, 

Page  122,  line  28: 

When  vertical  diffusive  transport  is  employed,  the 
computation  of  horizontal  parcel  advection  is  based  on  the 

Page  123,  insert  between  the  second  and  third  paragraphs: 

For  t  _ ns port  via  simple  advection  nlus  settling, 

TRAN?  is  cal  .  by  subroutine  ADVEC  with  XRIP=0  (otherwise 
KRIP=1) .  The  parcel  top  or  bittern  is  transported  stepwise 
through  the  vertical  layers  between  ZP  and  ZDEP  via  the 
layerwise  mode  (see  below).  In  each  layer,  with  index  KBH, 
the  vertical  velocity  is 
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WBAR  =  WFZ(KBH,  NDATO,  ..TIM)  -  CAVS(KBH) 


and  KAY  =  -  .1.  TSEG  is  computed  as  shown  by  the-  equation 
in  the  preceding  paragraph. 

Page  123,  third  paragraph,  line  3: 

is  traversed  at  a  time.  This  mode  is  mandatory  for  transport 
via  simple  adveetion  plus  settling,  or  when  a  parcel  trajectory 

Page  125,  line  6: 

The  rapid  computation  mode  is  employed  for  vertical 
diffusive  transport  when 

Page  145,  insert  at  the  end  of  the  Input  Data  Card  4  discussion: 

To  preempt  vertical  diffusive  transport,  set  KX=1  and 
set  ZMAX  arbitrarily  large.  This  causes  all  parcels  to  be 
transported  via  the  siraple-advection-plus-settling  mode 
(Eq.  (32)).  Of  course,  horizontal  diffusive  growth  of  parcels 
is  accounted  for  in  any  case. 

Page  146,  line  14: 

CSKIP  -  0.1 

Page  152,  line  8: 

HITIME.  If  the  vertical  diffusive  mode  of  transport  is  used, 
the  KBHX'th  base  should  be  above  or  at  the  top  of  the  trans¬ 
port  space  as  this  top  is  specified  by  ZMAX. 

Page  156,  Table  4,  Record  Number  12,  lir.e  $  under  Content: 
size  class  central  diameter  (gm),  mass  of  fallout  (kg) 

Page  154,  to  the  end  of  paragraph  1  add: 

However,  the  turbulent  energy  dissipation  rates  can  be  input 
only  for  the  horizontal  directions;  for  the  veicicel  direction 
Fickiar.  dif fusivities  always  are  input,  regardless  of  which 
typa  or  dac?  are  input  for  the  horizontal. 
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4.  CODE  REVISIONS 


4.1  Single  Card  Changes 

Addition  cf  the  arrays  CAVS(KSHF)  and  WAVG (KBhF,LTIMF) 
requires  revisions  in  DIMENS ION  statements  and  subroutine  argument 
lists.  However,  complete  FORTRAN  statement  listings  are  given  in 
this  supplement  for  all  subroutines  that  requite  these  revisions. 

Subroutine  EOUN,  card  27: 

CALL  NEST  (NET,  NETS1',  XO,  YO,  NDATO,  XL.  XR,  YL,  YU,  ICF, 
JCF,  NCF) 


Subroutine  DUMPER: 

Place  card  33  in  its  proper  position. 

Subroutine  NEST,  insert  between  cards  20  and  2'l 
DIMENSION  NET (ICF,  JCF),  NETSU(NCF) 


4.2  FORTRAN  Statement  Listings 

Complete  FORTRAN  statement  listings  are  given  for  the 
following  subroutines.  These  subroutines  are  operational  on  the 
UNIVriC  1108. 

Subroutine  Page 

C31M  i3 

ADMIN  18 

ADVEC  21 

AM3KT  23 

SPRVS  26 

TRANP  32 
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Ths  machine  used  to  prepare  these  listings  prints  a 
i ?  symto?  to  represent  a  4-8  punch;  this  symbol  should  be  an 
apostrophe  (’).  In  FORMAT  and  DATA  statements,  the  apostrophe  is 
used  to  define  Hollerith  character  fields. 


12 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

V 

c 

c 

c 

c 

c 

c 

“•I 

o 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


TRANSPORT 
IN  C31H. 


SEPTEMBER  1971 

C31M  IS  THE  MAIN  PROGRAM  WHICH  DIRECTS  THE  OIKFUSiyE 
MODULE  OPERATIONS. THE  OBJECT-TIME  DIMENSIONS  ARE  SET 
THESE  DIMENSIONS  ANO  THEIR  RESPECTIVE  ARRAYS  ARE 
«F  -  AA,BB,CC,DESOM,Evr  ,Q,H,OIFF,2HT 
LTIMF  -  TIHUP,OAVG,WAVG 
<8HF  -  Z8H,ZCH 
NOATF  -  KTOFO 

<BHF,  NDATF,  LTIMF  -  OFZ,OXSUN, OYSUM,USUM ,  VSl'M  ,  WFZ,  RSUH 

ICF,  JCF  -  NET 

NCR  -  NETSU 

MARF  -  MARY 

NAIF  -  ALT »ATEMP»RFO 

DATA  LITERALS  MUST  BE  INSERTED  IN  THE  DIMENSION  STATEMENTS  AND 
THE  RIGHT  HANO  SIDES  OF  THE  ARITHMETIC  STATEMENTS  IN  WHICH  THE 
ABOVE  VARIABLE 


IN 


AA(K) 

ALT 
ATEMP 

aeto 

CAV 

CAVS 
CC  ( K* 

CROSS 
CSKIP 
DA  VG 
OENCMOO  - 
DEP 

DIFFCK? 
OFKXS1- 
OF  l 

OTNCR  - 
DOEEN  - 
JOWM  - 
DT 

DWAF  - 
OXSUh  - 
SYSUH  - 
CZ 

DZMIN  - 

E<K*  - 
EODY  - 

EFFLUX 
F{<*  - 

FAV 

FMAB  - 
FMBEL  - 
TCF 
ICX 

TMARIN- 
TPOL’T  - 
ISIN  - 


NAMES  APPEAR. 

GLOSSARY  *»*****»**»***7**»**»U***** 

EQUALS  S2*  (CTFFfK*l»+OIFFmi  +  Sl*F{K«-l* 

ALTITUDES  FCR  ATMOS.  DENSITY  AND  VISCOSITY  TABLE 
DYNAMIC  VISCOSITY  OF  AIR  DATA  VECTOP  FOR  ATMOS.  TABLE 
EQUALS  S2MCIFFf<*ii+2.*DIFFj<)4-0IFFf<-ll)«-Sl*F{<) 

AVG.  FALLRATE  USED  IN  COMPUTING  POEST .  IT  APPLIES  MID  WAY 
FROM  PARCEL  TOC  TO  ZMIN. 

PARTICLE  FALL  RATE  FOR  EACH  ATMOS.  STRATUM 
EQUALS  SIMCIFF«K)*-OIFFIK-m 

CROSSWIND  CROSSING  TRAJECTORIES  CCRRECTION  TO  TURB. 

TOTAL  FRACTIONAL  PARCEL  DEPOSITION  THRESHOLD 
AVG.  ATMOS.  VERT.  TURB.  PER  UPDATE  DATA  VECTOR 
EQUALS  1.  ■  THETA *(BB(K*-CC  fKI*E  CK-1) * 

DEPOSITED  FRACTIONAL  MASS  INCREMENT 

-  VERT.  OIFFUSIVITY  AT  K-TH  SHALL  ALTITUDE  INCREMENT 
VERTICAL  OFFUSIVITY  AT  ALTITUDE  INCREMENT  <k-l 
TURBULENCE  Z  COMPONENT  3-OIM.  DATA  ARRAY 
RATE  OF  CHANGE  OF  FRACTIONAL  MASS  DEPOS.  RATE  THRESH. 

MASS  DEPOSITION  RATE  THRESHOLD 

DOWNWIND  CROSSING  TRAJECTORIES  CORRECTION  TO  TURB. 

SMALL  ITERATION  TINE  STEP  FOR  VERT.  OTFF.  DIFF.  EQ. 

PARCEL  VERT.  THICKNESS  BEFORE  ADVECTION 
TURBULENCE  V  COMPONENT  CWEIGHTEO  SUM!  3-DIM.  DATA  ARRAY 
TURBULENCE  Y  COMPONENT  (WEIGHTED  SUM*  3-DIM.  DATA  ARRAY 
SMALL  ALTITLOE  INCREMENT  FOR  VERT.  DIFF.  OIFF.  EQ. 

MINIMUM  VALLE  OF  DZ 
EQUALS  THETA*AA1KI/OENOHIK» 

RATIO  OF  LAGRANGIAN  TURBULENCE  TIME  SCALE  TO  EULERIAN 
TUR8ULENCE  LENGTH  SCALE 
UPPER  EFFLUX  FRACTIONAL  MASS 

IN  SUB.  DIFFF,  WORKING  SPACE  IMPLICIT  METHOD  DATA  VECTOR. 
IN  SUB.  AM3NT,  WORKING  SPACE  pOR  VERTICAL  VELOCITIES 
MID-ATMOS.  AVG.  FALLRATE.  USED  IN  CROSSING  TRAJECTORIES 
CORRECTIONS  AND  IN  TRUNCATION  ERROR  ESTIMATION  . 
CUMULATIVE  FRACTIONAL  MASS  AIRBORNE 
MIN.  PARCEL  FRACTIONAL  MASS  ALOFT  TO  BE  TRANSPORTED 
MAX.  FORMAL  DIM.  CORRESPONDING  TO  ICX 

OOJECT-TIHE  FIRST  MAX.  DIM.  OF  ARRAY  NET.  NUMBER  OF  NET 
MESHES  IN  EAST-WEST  ROW. 

LOGICAL  UNIT  NUMBER  OP  CR-TRSNS.  INTER.  MOD.  OUTPUT  TAPE 
LOGICAL  UNIT  NUMBER  OF  CIFF.  TRANS.  MOD.  OUTPUT  TAPE 
LOGICAL  UNIT  NUMBER  OF  SYSTEM  INPUT  TAPE 


C31M 
C31M 
C31M 
C31M 
C31H 
C3  j  M 
C31M 
C31M 
C31M 
C31M 
C31M 
03 1M 
C31M 
C31M 
C31M 
C31H 
>C31H 
C31M 
C31M 
C31M 

C31M 

C31H 

C31M 

C31M 

C31M 

C31M 

C31M 

C31M 

C31M 

C31M 

C31M 

C31H 

C31M 

C31M 

C31M 

C31M 

C31M 

C31M 

C31M 

C31M 

C31M 

C3LM 

C31H 

C31M 

C31M 

C31M 

C31M 

C31M 

C31M 

C31M 

C31M 

C31M 

C3IM 

C31M  - 

C31H 

C31M 

C31M 

C?1  M 


1 

? 

3 

4 

5 

6 

7 

8 
B 

10 

11 

12 

13 

14 

15 
If 

17 

18 
IP 
20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 
3? 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 
5% 

55 

56 

57 

58 
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ooooaooooooooocoooooonoooooooooooooooooooooooooooooooooooo 


532^, 


isour 

JCF 

JCX 


<8H 


<BHF 

K8HX 

«F 

KKM 

KKMfti 

«X 

KKXS1 

<«>IP 


KTOPO 

KX 

<XMIN 

LSTEP 


LTIM 

LTTMX 

LTIMF 

MARF 

MARX 

MA°Y 


HC 


LOGICAL  UNIT  NUMBFR  OF  SYSTEM  OUTPUT  TAPE 
MAX.  FORMAL  DIM.  CORRESPONDING  TO  JCX 

OBJECT-TIME  SECOM)  MAX.  DIM.  OF  ARRAY  NET.  NUMBER  OF  NET 
MESHES  IN  SCUTH-NORTH  ROM. 

ATMOS.  VERT.  SPACE  INDEX  FOR  ARRAYS  USUH,VSUM, WFZ,DXSUH , 

OYSUM,  OFZ,  RSUH,  ZBH,  ZCH 

MAX.  FORMAL  DIM.  CORRESPONDING  TO  <BHX 

03JECT-TIME  MAX.  VALUE  OF  KBH 

MAX.  FORMAL  OIM.  CORRESPONDING  TO  «X 

ALWAYS  EQUALS  2.  CORRESPONDS  *G  K=0  ALTITUDE  INCREMENT 
EQUALS  KKMM 
EQUALS  <X«-<XM 
EQUALS  KKX- 1 
•  CONTROL  VARIABLE 

0  FOR  AO VEC T I VE  TRANSPORT 
1  FOR  DIFFUSIVE  TRANSPORT 

NET  MESH  ANC  SUP-MESH  TOPOGRAPHY  TABLE  DATA  VECTOR 

MAX.  NUMBER  OF  OZ  ALTITUDE  INCREMENTS 

MIN.  NUMBER  OF  DZ  ALTITUDE  INCREMENTS 

NUMBER  OF  IMPLICIT  METHOD  ITERATIONS.  n-eLT=LSTEs*OT. 

DYSUM,  OFZ,  RSUH,  TIMUP,  0 AVG,  WAVG 

ATMOS.  UPDATE  TIME  INDEX  FOR  ARRAYS  ”SUM,VSUM»WFZ,DXSUM, 

OBJECT-TIME  MAX.  VALUE  OF  LTIM 

MAX.  FORMAL  OIM.  CORRESPONDING  TO  LTIMX 

MAX.  FORMAL  OIM.  CORRESPONDING  TO  MARX 

OBJECT-TIKE  MAX.  DIM.  OF  ARRAY  MARY 

HORIZ.  ATHCS.  SPACE  RESOLUTION  NET  MESH  AND  SUB-MESH 
CONTROL  FLAGS  DATA  VECTOR 


CONTROL  INTEGER  DATA  VECTOR 


KCC1)  LESS  THAN  OR  EQUAL  TO  ZERO,  SUPPRESSES  LISTING  OF 
EXPANOED  WIND  AND  TURB.  DATA 


MC ( II  GREATER  THAN  OR  EQUAL  TO  ONE,  CAUSFS  LISTING  OF 
WINO  AND  TURB.  OATA  BEFORE  SUMMATION 


MC ( II  EQUALS  TWO,  CAUSES  LISTING  CF  KIND  AND  TURB.  DATA 
AFTER  SUMMATION 


HCC21  EQUALS  ONE,  SUPPRESSES  LISTING  OF  ATMOS.  VISC.  AND 
DENS.  TABLES 


MC (3)  EQUALS  ZERO,  SUPPRESSES  LISTING  OF  OEPOSTT 
INCREMENTS  CN  TAPE  ISOUT 


HC  f 4)  EQUALS  ONE,  CAUSFS  PRINTOUT  OF  TRANSPORT 
INTERMEDIATE  RESULTS  ON  TAPE  ISOUT.  WARNING.  PRINTOUT 
IS  EXTRAORDINARILY  VOLUMINOUS.  FOR  DEBUGGING  ONLY. 


HCC7)  EQUALS  ONE,  SUPPRESSES  LISTING  OF  RAW  WINO  AND  TURB. 
INPUT  DATA 


NCC10I  EQUALS  ONE,  CAUSES  TURB.  DATA  TO  BE  TREATED  AS 
KOLMOGOROFF-BATCHELOR  ENERGY  DISSIPATION  RATES 


MC  CIO I  NOT  EQUAL  TO  ONE,  CAUSES  TURB.  DATA  TO  BE  TREATED 
AS  FICKIAN  OIFFUSIVITISS 


C31M 

59 

C31K 

60 

03 1H 

61 

C31M 

62 

C31M 

63 

C31M 

64 

C31K 

65 

C31M 

66 

C31M 

6? 

C31M 

66 

C31M 

69 

C31M 

70 

C31M 

71 

C31« 

?2 

C31M 

7? 

C31M 

74 

C31M 

75 

03 1M 

76 

C31M 

77 

C31M 

78 

03 1 M 

fo 

C31M 

80 

C31M 

61 

C31M 

«2 

C31M 

53 

C31M 

8% 

031M 

85 

C31M 

86 

C31H 

87 

C31M 

55 

C31M 

89 

C31M 

90 

C31M 

91 

C31K 

92 

C31M 

9? 

C31« 

94 

031 H 

95 

C31M 

96 

C31M 

97 

C31M 

96 

C31M 

99 

C31M 

100 

C31M 

101 

C31M 

102 

C31M 

103 

C31M 

104 

C31M 

1C5 

C31M 

106 

C31M 

107 

C31M 

108 

C31M 

109 

C31M 

110 

C31M 

111 

C31M 

112 

C3iM 

113 

C31M 

114 

C31M 

115 

C31M 

116 
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oooooooooooc: 


MINT  - 
NAT 

NATE  - 
NRL<  - 
NCE 
NCX 

NDATA  - 


NOATC 

NDATO 

NOATP 

NOATX 

NOATE 

NDELT 

NET 

NETSU 

NSEQC 

PHI 

3(K) 

RHO 

PC 

ROPART 

RSUH 

RWAF 

SIGXG 

SIGYC 

51 

52 

TOELT 

TOEP 

THFTA 

THETQ 

TIKE 

TIMEX 

TIMUP 

TO 

TP 

TPAUS 

USUM 

VETA 

VSUM 

WHO 

WAVG 

HAVGK 

WF? 

HINT 

XLLC 

xo 

XP 

YLLC 

YO 

YP 

ZfiH 


MC C 18)  EQUALS  CNF,  SUPPRESSES  READING  FROM  TAPE  IPARIN 
ANC  WRITING  CNTC  TAPE  IPOUT 


MIN.  HUMBER  OF  CT  SMALL  TIME  STEPS  PER  CEPOSIT  INCRFHENT 
TIME  INTERVAL 

NUMBER  OF  ALTI T'JOE  STRATA  IN  ATMOS,  OENS.  AND  VISC.  TABLE 

MAX 3  FORMAL  DIM.  CORRESPONDING  TO  NATF 

RECORD  BLOCK  SIZE  FOR  DEPOSIT  INCREMENT  »EC0»03 

MAX .  EORMAL  DIM.  CORRESPONDING  TO  NCX 

OBJECT  TTNE  MAX.  DIM.  OF  AtTRAY  NETSU 

AfPU?.  HORI ?«  SPACE  INDEX  FCP  ARRAYS  USUMt VSUM,WEZ, OXSUM , 
OYSUN,  OFZ,  RSUM,  <TOPO 

HORIZONTAL  INDEX  OF  LATTICE  CELL  CONTAINING  POINT  (XCtYC) 

HORIZONfAL  INDEX  OF  LATTICE  CELL  CONTAINING  POINT  CXO,YO) 

HORIZONTAL  INDEX  OF  LATTICE  CELL  CONTAINING  POINT  (XP,Y») 

OBJECT -T IME  MAX.  VALUE  OF  NO AT  A 

MAX.  FORMAL  DIM.  CORRESPONDING  TO  ND4TY 

NOMINAL  NUMBER  OF  DEPOSIT  INCREMENTS  PER  FALLOUT  PAPCEL 


VISC.  TABLE 


HORIZON!  AL  INDEX  t 
HORIZONTAL  INDEX  ! 
OBJECT-TIME  MAX.  \ 
MAX.  FORMAL  DIM.  f 
NOMINAL  NUMBER  OF 


C31* 

C31M 
C31 M 
C3 1 M 
C31M 
C31M 
CS1* 
C3 1M 
C31M 
C31M 
C3 1M 
C31M 
C31M 
C31M 
C31M 
C3i* 
C3  1M 
C31M 
C31M 
C31  M 
C3H 


HORIZONTAL  SPACE  CONTROL  NET  MESH  2-DTM.  ARRAY  C31M 
HORIZONTAL  SPACE  CONTROL  NET  SUB-MESH  DATA  VECTOR  C3H 
STORAGE  SEQUENCE  ORDINAL  OF  FIRST  PARCEL  TO  BF  TRANSPGRTF0C31M 
EQUALS  1-THETA  C31M 
CONCENTRAT  T CN  IN  K-TH  ALTITUDE  INCREMENT  C31M 
ATriOS.  DENSITY  DATA  VECTOR  FOR  ATMOS.  TABLE  C3JM 
KINO  HEADING  ORIENTATION  ANGLE  AFTER  ADVECTION  C3SH 
FALLOUT  PARTICLE  DENSITY  C31M 
HIND  HEADING  ORIENTATION  ANGLE  (WEIGNTED  SUM)  3-OIM.  APR AYC31M 


PARCEL  RAOILS  IN  PARCEL  CENTRAL  PLANE  8EF0RE  ADVECTTON 

FARCEL  MASS  H0«o  STAND.  DEV.  DOWNWIND  AFTER  ADVECTION 

PARCEL  MASS  HOR,  STAND.  DEV.  CROSSKTND  AFTER  ADVECTION 

EQUALS  OT/OZ 

EQUALS  OT/ C2.*(DZ)**2)I 

CURRENT  OEPCiXT  INCREMENT  TIME  INTERVAL 

ADVECTIVE  T FAN SPORT  TIME  INTERVAL 

IMPLICIT  FIMTE  "INFERENCE  PARAMETER 

OOU9LE  PRECISION  WORD  CORRESPONDING  TO  THETA 

TIME  AT  ONSFT  CF  CURRENT  DEPOSIT  INCREMENT  TIME  I NTCRV AL 

SIMULATED  TRANSPORT  TIME  LIMIT 

ATMOSPHERE  LPDATF  TIMETABLE  DATA  VECTOR 

TIME  AFTER  PARCEL  A0VECTI3N 

TIME  BEFORE  PARCEL  ADVECTION 

TIME  AT  ENG  CF  CURRENT  DEPOSIT  INCREMENT  TIME  INTERVAL 

WIND  X  COMPONENT  (WEIGHTED  SUM)  3-DIM.  DATA  ARRAY 

VERTICAL  DIFFUSION  ABSORPTION  COEFFICIENT 

HIND  Y  COMPONENT  (WEIGHTED  SUM)  3-OIM.  DATA  ARRAY 

SETTLING  RATE  AT  K-TH  SMALL  ALTITUDE  INCREMENT 

AVG.  ATMOS.  VE°T .  WIND  PER  UPDATE  PER  STRATUM 

KAVG  AVERAGED  OVER  THE  STRATA  FOR  THE  FIRST  UPDATE 

HIND  Z  COMPONENT  3-THH*  DATA  ARRAY 

NET  CONTROL  MESH  DIMENSION 

X  COORDINATE  OF  SOUTH-WEST  CORNER  OK  ATMOS.  SPACE 
PARCEL  CENTER  X  COORDINATE  AFTER  ADVECTION 
PARCEL  CENTER  X  COORDINATE  BEFORE  AOVEGf ION 
X  YOOROINATE  OF  SOUTH-WEST  CORNER  OF  ATMOS.  SPACE 
PARCEL  CENTER  Y  COORDINATE  AFTER  ADVECTION 
PAPCEL  CENTER  Y  COORDINATE  BEFORE  ADVECTION 


COORDINATE 

COOROINATi- 


-  ATMOSPHERE  STRATA  BASE-ALf HUf'E  DATA  VECTOR 


C31M 
C31M 
C31M 
C31M 
C31M 
C31M 
C31M 
C31M 
C31M 
C31M 
C31M 
C31M 
C31M 
C31M 
C31M 
C3 1M 
C31M 
C31M 
C31M 
C31« 
C31M 
C31M 
C31M 
C31H 
Coin 
C31M 
r.3i« 
C31H 
C31M 
C3  IS 
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ZCH  -  ATMOSPHERE  STRATA  MID-ALTITUDE  DATA  VECTOR 
ZOEP  -  AOVECTIVE  TRANSPORT  TERMINAL  ALTITUDE 
ZHT(K»  -  K-TH  ALTITUDE  INCREMENT  ABOVE  7MIN  LEVEL 
ZLOW  -  PARCEL  BASE  ALTITUOE  BEFORE  ADVECTION 

ZMAX  -  ATMOSPHERE  TOP  ALTITUDE  IFOR  VERT.  OIFF.  DIFF.  EQ.3 

7MIN  -  GROUND  LEVEL  OR  DEPOSITION  PLANE  ALTITUDE 

ZO  -  PARCEL  CENTER  Z  COORDINATE  AFTER  ADVECTION 

ZP  -  PARCEL  CENTER  Z  COORDINATE  BEFORE  ADVECTION,  ECXEPT  AS 

REDEFINED  IK  SUP.  ADVEC 

ZUPP  -  PARCEL  TOP  ALTITUDE  BEFORE  ADVECTION.  ZLOW+CWAF 

DOUBLE  PRECISION  A  A  C  2  0  *t  5  ,8«t  2043  ,CC  (204)  ,  DENQM  (2  041  ,E C704 ) ,F f 2 04 > 

DOUBLE  PRECISION  01204) 

DIMENSION  Z9H<  27),ZCH<  27),TIMUP<  6),MARY(  i» 

DIMENSION  KTOPOC  t),NETSU<  1) ,  NET !  i,  13 


DIMENSION  Z9H< 
UIMENSiON  KTOPOi 
DIMENSION  DFZC 
DIMENSION  USUM { 
DIMENSION  OXSUMf 
DIMENSION  RSUH { 
DIMENSION  WAVGi 


6) ,  WFZ i 
6) ,  VSUH( 
63 , D YSUM  C 
6) ,0AVG< 


6) , MARY ( 
1* 

27, 

27 , 
27, 

6) 


NATF=260 
NCF=i 
NO ATF=1 
DO  1  N=1,NCF 

1  NETSU (N) =0 
DO  2  J=1 , JCF 
DO  2  1=1, ICF 

2  NET(I,J1=0 

DO  3  M=1 ,M ARF 

1  AAA  rtU  lUl.n 

y*  t  I  r%  ‘\  1  %  I  I  f  —  U 

00  4  <=1;K8HF 
ZputMisO. 

4  ZCHf<i^S. 

D'"-  10  6  S=i.K5K? 
DO  104  L=i ,lTIhF 
iOE  «A¥GtiC,D-3.G 
DO  5  n=1,NOATF 

5  KTOPO (N3  =0 

00  6  L=1,LTIMF 

TIMUPa»=0. 

DAVG(L>=0. 


C31M 

C31M 

C31U 

C31M 

C31M 

C31M 

C31M 

C31M 

C31H 

C31M 

C31M 

C31« 

C31M 

C31M 

C31M 

CUM 

C31M 

CM* 

C31M 


DIMENSION  ALT{  260) 

,RHOt  260), 

ATEMP ( 

2603 

C31  * 

194 

DIMENSION  CAVSl  27» 

,Wt  204) ,DIFF  C  204) »ZH?  { 

204) 

C31M 

199 

COMMON  /QPARM/  IPOUT 

,IPARIN,NBLK 

,NAT 

»NDELT 

,*X 

,  K KM 

C31M 

196 

1 , NSEQO  ,ICX  ,JCX 

, NCT  , KBHX 

vNDATX 

,LTIMX 

,  ISIN 

, ISOUT 

C31M 

197 

2, EDDY  , FMREL  ,L$?EF 

,NCtl3) , HINT 

,  XLLC 

,  YLL  C 

,yhfta 

•  ZMI N 

C31M 

193 

3,CS<IP  , MINT  , ZMA X 

, TIMEX  ,OT 

♦  DZ 

,X° 

,  YP 

»2P 

C31M 

199 

4, 01 NCR  , DOWN  , TP 

,  ZLOW  , OWAF 

,  PNAF 

, POPART 

,ZU  PP 

,  VET  A 

C31M 

200 

5 , DOpEN  , CROSS  , T tME 

,<i<MAl  ,«X 

,«<S1 

,<XMIN 

,NOATP 

CMM 

201 

ICF=  1 

C31M 

202 

ISIN=5 

C31M 

203 

ISOUT=6 

C31M 

204 

IPARIN=9 

C31M 

209 

IPOUT=10 

”31M 

206 

JCF  =  1 

C31M 

207 

<BHF=27 

C31M 

20  B 

KKF=204 

C31M 

209 

LTIMF=6 

C31M 

210 

MARF=i 

C31M 

211 

C3 1M 
C31M 
C3  1M 
C31M 
C31M 
C31M 
C31M 
C31M 
C31M 
C31M 
C31M 
C31M 
C31M 
C31M 
C31M 
C31K 
C31M 
C31M 
C31M 
C31M 
C31M 


on  6  N  =  1  »NC ATF  C?1M 

DO  6  rf=l,K9HF  C31* 

DFZ  (K,N,L»  =  0.  C3  }M 

WFZ (<  »N»L) =0.  C3JM 

USUH<K»N»LJ-0«  03 1M 

VSUMC<,N,U=0.  C3  JH 

OXSUM(<,N,U  =  0»  C31M 

DYSUM(<,N,L>=0.  03 1M 

6  RSUH<K»N,L>=0.  031M 

COMMENCE  READING  DATA  INPUTS  FROM  TAPES  ISIN  AND  IPARIN  C31« 

COMMENCE  WRITING  OATA  OUTPUT  HEADERS  ONTO  TAP^S  ISOUT  AND  IPOUT  C31M 

CALL  LIN<{ALT  ,RHO?AT£MP,NATF)  C31M 

CONSTRUCT  THE  HORIZONTAL  SPACE  CONTROL  NET  C31« 

CALL  GETU»*(NET,NETSU,KTOPC,HASYfMARF,ICF,  JCF.NCF, NOATF)  C31M 

CONSTRUCT  THE  ATMOSPHERIC  STRATA  ANO  UPDATE  DATA  VECTORS  CM* 

CALL  HITIHF(ZCU,ZOH,TIHUP,KPHF,LTIHF»  C31M 

CONTSRUCT  AND  FILL  IN  THE  ATMOSPHERIC  LATTICE  AND  UPDATE  STRUCTURE  C3 1 H 

CALL  AnMIN(NET,NETSU,ZnH,ZCH, TIHUP,USUH, VSUU,nXSUH,OYSUH,  C31M 

1RSUM,0FZ,WFZ,0AVG, WAVG, ICC, JCF, NCF , X8HF ,NOATF ,LT IMF)  C31H 

CIRCUMVENT  ALL  TAPE  HANDLING  IF  MCC18*  EQUALS  1  C31M 


IF«MCC13).ea.ll  GO  TO  7 

CALCULATE  THE  OIFFUSIVF  TRANSPORT  OF  PARCELS  ACCEPTED  FROM  TAPE  IPAPIN 
COPY  CUT  RESULTS  ONTO  TAPE  IPOUT 

CALL  SPRVS(NET,NETS'J,Z8H,ZCH,TIMUPtUSUH,  VSUM,DXSUM,OYSUM, 

1RSUM,DFZ»WFZ»DAVG»NAVG,ALT»RH0»AT£NP,AA,BB»CC»DEN0M,PIFF,£,F,Q,W  , 
2ZHT,  ICF,JCF,NCF,(<BHFtNCATF,LTlMF,(<<F,NATF,CAVS» 

7  CALL  EXIT 
STOP 
END 


C31M  25? 
C31M  254 
C31H  255 
C3 1M  256 
C31M  257 
C31M  25? 
C31M  259 
C31H  260 
C31M  261 


SUBROUTINE  AOMIN (NET ,NETSU»ZBH,ZCH»TIMUP,USUM,  VSUM,DXSUH, DYSUM, 
1RSUM, OPZ,WFZ,DAVG,W AVG, ICF,JCF,NCF, KBHF, NOATF,LT IMF) 

SEPTEMBER  19/1 

SUBROUTINE  AOMIN  CONSTRUCTS  NINO  DATA  ARRAYS 
USUM,  VSLM,  HFZ,  WAVG,  RSUH 
ANO  TURBULENCE  DATA  ARRAYS 

QXSUM,  DYSUM,  OFZ,  OAVG . 

IN  AOMIN  ONLY  LTIM  ANO  SPEC  ARE  RE  AO  FROM  TARE  ISIN. 

SPEC  -  OATA  SPECIES  IDENTIFICATION  WORO  *WTNO*  OP  *SPEC* 

LTIM  -  UPDATE  ORDINAL  OF  OATA  SET.  FIRST  ATMOS.  SET  HAS  LTIH=i, 
SUB.  HKOAT  IS  CALLEC  TO  PERFORM  DATA  EXTRAPOLATIONS. 

AREA  -  AREA  OF  HORIZ.  SPACE  NET 

AREAN  -  AREA  OF  N-TP  NET  MESH  OR  SUB-MESH 


W 


ADMIN  1 
AOMIN  ? 
ADMIN  3 
AOMIN  4 
AOMIN  5 
AOMIN  6 
AnpiN  7 
AOMIN  F 
AOMIN  9 
SDMTN  10 
AOMIN  11 
AOMIN  12 
ADMIN  13 
AOMIN  14 


COMMON 

/QPARH/ 

IPOUT 

,  IPAR IN , NBLK 

,NAT 

, NDELT 

,<X 

,«M 

AOMIN 

15 

1.NSEQ0 

,ICX 

,  JCX 

, NCX  , KBHX 

,NOATX 

,LT IMX 

,  ISIN 

, ISOUT 

ADMIN 

16 

2 , EDDY 

, FMBFL 

,LSTEP 

, MC(18) , WINT 

,  XLLC 

,  YLLC 

, THETA 

,  ZMTN 

ADMIN 

17 

3,CS<IP 

,  MINT 

,  ZMA  X 

, TIMEX  ,DT 

» 07 

,XP 

»YP 

,7P 

ADMIN 

Id 

4, DINT  R 

,OOWiJ 

,TP 

, ZLDW  , DWAF 

,  RWAF 

,ROpART 

,ZUPP 

,  VET  A 

AOMIN 

19 

5,0 OPEN 

,  CROSS 

,  TIME 

, KKMA 1  , KKX 

,<KXS1 

,KXMIN 

,NOATP 

AOMIN 

20 

DIMENSION  RSUM{KBHF,NOATF,LTIMFI,OAVGCLTIMF)9WAVG(KBHF,LTIMF»  AOMIN  21 

DIMENSION  NET ( ICF, JCF) ,NETSU(NCF) ,7CH(<OHF) ,TIMUPfLTlMF) ,ZBH(<9HF>  AOMIN  22 
OIMENSION  USUM (KBHF  ,NO ATF,LT IMF) , V  SUM ( KBHF , NQA  TF , LTIMF)  ADMIN  23 

DIMENSION  OXSUMf  K8HF,  NDATF  ,LT  IMF) , OYSUM(KBHF,NDATF,LT IMF)  AOMIN  24 

DIMENSION  OFZ (KBHF, FOATF,LTIHF) ,WFZ (KBHF, NOATF ,LT IMF)  AOMIN  25 

DIMENSION  LH(10) ,LO  (10)  AOMIN  26 

INTEGER  WlNO»QFSN, 0CNE,SpEC  AOMIN  27 

OATA  PROGRM/*ADMIN  */  AOMIN  23 

DATA  WlND/*WlND*/  AOMIN  29 

OATA  OFSNY*OIFF*/  AOMIN  30 

OATA  DONE/*NO  M*/  ADMIN  31 

FORMAT(*0*36X, Y UPDATE  TIME  INDEX*I5»*.  WlNO  GRIO  CELL  IN0EX*I5/)  AOMIN  32 
FORMAT ( *  WIND*2(6X,*HCRIZCNTAL*) ,6X,*  VERTICAL  *6X,*CROSSWINO  *6XADMIN  33 
1, *  OOWNWINO  *6X,*  VERTICAL  *6X,  HORIZONTAL*)  AnMTN  34 

3  FORMAT ( *  LAVER*6X, <E.-H.  WTN0*6X,*N.~S.  WINDE6X, *  HIND  *3(6X,*ADMIN  35 

20IFFUSI0N  *),6X,*  RCTATION*)  AOMIN  3F 

4  FORMAT (*  IN0EX*3(6X,*  VFLOCITY  *),3(6X,*  CONSTANT  *),6X,*  ANGLE*) AOMIN  37 

5  FORMAT <*  *I5,7E16.4)  AOMIN  33 

6  FORMAT (/ 25 X,* WEIGHTED  SUMS  OVER  ABOVE  COLUMN  ENTRIES*/)  ADMIN  39 

7  FORMAT f *  *I5,2E16.4,16X,2E16.4,16X,E16.4)  AOMIN  40 

8  FORMAT (*  *22X,*UPDATE*I4,*»  CELL*I4,  **  AOMIN  41 

3AVG.  VERT.  OIFF.  =*E12.4)  AOMIN  42 

9  FORMAT (56X,I2,8X,3A4)  ADMIN  43 

10  FORMAT (*25X,*ATM0SPFERE  UP0ATE*I4,*  FOR  TIMES  LATER  THAN  *E12.4,*  AOMIN  44 

1SECONOS*)  AOMIN  45 

11  FORMAT  It  *25X ,  **  *********  WINDFIELO  EXTRAPOLATION  *  *  *  *  AOMIN  46 

1  *  *  *  *  **/)  AOMIN  47 

12  FORMAT (*  *25X,**  '*******  OIFFUSIVITY  EXTRAPOLATION  *  *  *  *AOMIN  48 

2  *  *  *  *  **/)  ADMIN  49 

13  FORMAT (/ 25 X,*UpO ATE *14,*  OF  THE  HINDFIELO  IS  MISSING*)  AOMIN  50 

14  FORMAT (/25X,*UP0ATE*I4,*  OF  THE  DISFUSIVITY  IS  MISSING*)  ADMIN  51 

15  FORMAT (*  OVER  ENTIRE  H0RI70NT AL  GRID  FOR  UP0ATE*T4,  AOMIN  52 

6  **  AVG.  VERT.  OIFF,  =*E12.4)  ADMIN  53 

16  FORMAT (*  AVG.  VERT.  VEL.  FOR  EACH  STRATUM  IS  -  *1  AOMIN  54 

17  FORMAT (2X, 15, E 16.41  AOMIN  55 

AREA  =  ICX*JCXMHTNT**2I  AOMIN  56 

00  999  L=1,LTIMX  AOMIN  57 

LW(L)=L  ADMIN  58 


OIFFUSIVITY  EXTRAPOLATION 
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999  LOCL) =L  ADMIN  59 

COPY  IN  cTIM  AND  SPEC  FROM  TAPF  ISIN  AND  CALL  SUB.  MKDAT  AOMIN  60 

1000  REA0;iSIN,9>  LTIM,SPEC  ADMTN  61 

IF (SPECeEQ.DCNEI  60  TO  1301  AOMIN  6 2 

IF{ (LTIH.LT.l) .OR. (LTIH.GT.LTIMXj )  CALL  ERROR t »POGRM,-100 0, TSOUT )  AOMIN  63 
WRITE(ISOUT,10)  LTIM,TIMUP!LTIM;  AOMIN  64 

IF ( SPEC. EQ. HINDI  GO  TO  1101  ADMIN  65 

IF(SPEC.EQ.OFSN)  GO  TO  1201  ADMIN  66 

CALL  ERROR (PROGRM,- 1101, ISOUT!  ADMIN  67 

1101  WRITE ( ISOUT  ,11 )  AOMIN  68 

00  1102  t=i,LTIMX  ADMIN  69 

1102  IF(LTIM.EQ.LH(LI)  GC  TO  1103  AOMIN  70 

CALL  ERROR  (PROGRM, -1102, ISOUT)  AOHIN  71 

1103  LWCL>=-1  ADMTN  7? 

CALL  HKDAT(ZC»-,NET,NETSU,LTIM,  USUM,  VSUM, WFZ, ICFt JCF, NCF,  AOMIN  72 

1K3HF,N0ATF,LTIHF)  ADMIN  74 

GO  TO  1000  ADMTN  75 

1201  WRITE (ISOUT  ,12)  AOMIN  76 

DO  1202  L-i ,LT IMX  ADMIN  77 

1202  IF(LTIrt.EQ.LO(Lll  GC  TO  1203  AOMIN  78 

CALL  ERRCMPROGRM, -1202, ISOUT)  ADMIN  79 

1203  LD (!. )  =-l  ADMIN  80 

CALL  MKOAT(ZC»-,NET,NETSU,LTIM,nxSUH,OYSUM,OF?,ICF, JCF, NCF,  ADMIN  51 

! <8HF» NOATF, LTIMF)  ADMIN  82 

GO  TO  1000  ADMIN  83 

CHECK  IF  ANY  WINO  DATA  SETS  APE  MISSING  ADMIN  84 

1301  00  1303  L=1,LTIMX  ADMIN  85 

IF(LWCL).EQ.-l)  GO  TO  1302  ADMTN  86 

WRITE(ISGUT,13I  LTIM  ADMIN  87 

CALL  ERROR (PROGRM, - 1302,ISOUT1  AOMIN  93 

CHECK  IF  ANY  TURBULENCE  DATA  SETS  ARE  MISSING  AOMIN  89 

1302  IF(LO(L).EQ.-l)  GO  TO  1303  ADMTN  90 

WRITE (ISOUT, 141  LTIM  ADMIN  Cj 

CALL  ERPORCPROGRM, -1303, ISOUT!  ADMIN  92 

1303  CONTINUE  ADMIN  93 

CALCULATE  THE  WEIGHTED  SUMS  OVER  ATMOS.  STRATA  AND  REWRITE  ARRAYS  ADMTN  94 

C  USUM,  VSUM,  RSUH,  OXSYM,  OYSUH.  ALSO  COMPUTE  DAVG  AND  WAVG.  ADMIN  95 

ZSPAN=ZBH(K6HX)-ZBH  (1)  ADMIN  96 

DO  922  L=1 ,LTIMX  AOMIN  97 

CO  1304  LK=1,K8HX  ADMIN  98 

1304  WAVG(LK,ll  =  Q.  0  ADMIN  99 

DAVG(L)=0.  ADMIN100 

DO  921  N=1 , NDA TX  ADMINlfJl 

DKAV=0.  ADMIN102 

IF  (MC  (1) .LT.l)  GO  TC  915  aDMIN103 

WRITE ( ISOUT ,11  L,N  ADMIN104 

WRITE ( IS OUT,2)  ADMIN105 

WRITE (ISOUT ,3)  ADMIN106 

WRITE (ISOUT, 41  ADMIN107 

915  DO  920  K=1 , KBHX  ADMTN10 9 

UKNL=USUM(K,N,U  ADMIN109 

VKNL-VSUMf K,N, Lt  ADMIN110 

IF(ABS(U<ND-1.0E-30)  9151,9154,9154  ADMIK'lll 

9151  IF (A8S (VKNLI-1.  OE-3  0!  9152,9153,9153  A0MIN112 

9152  R<NL=0.  AOMIN113 

GO  TO  9155  ACMIN114 

9153  RKNL=1, 57079633  A0MTN115 

GO  TO  9155  A0MIN116 
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9154  RKNL=ATAN(VKNL/U<NL1 

9155  DXKNL=OXSUM(K,N,L) 

OYKNL=OYSUM  1K,N,LI 

IF (<-<8HX>  916,9155,920 

9156  R5UM«K,N,L»=RKNL 
GO  TO  9195 

916  ZSTEP=ZBH(K«-1>-Z9H(X> 

USUMCK,N,L*=UXNL*ZSTEP 

VSUM(<,N*LI=V<NL*Z5TEP 

RSUM(K,N,LI=RKNL*ZSTEP 

0X?UM(<,N,LI=9X<NL*ZSTEP 

OYS!JM(K,N,L)  =  OYKNL*ZSTFP 

OSU^-ZST£P*OFZ(<,N,U 

M=K’-1 

IF (Ml  920,9i9, 918 

918  USUM(K,N,L)=U9UM(^,  K,L)  «-USUM(M,N,L» 
V5UMU,N,L)=VSUM(<,h,L)*-VSUMtH,N,L) 
RSUM(K,N,L)=RSU«(K,*,L)4-RSUM{M,N,L) 

OXSUH (K ,N»  LI=OXSUM (K,N,L)+OXSUH(H,N,L) 
DYSUM(<,N,L>=OYSUH(K,N,L)+OYSUM(M,N,L> 

919  OKAV=DKAV*DSUM 

9195  IF(MC(1) .LT.l)  GO  TC  920 

WRITE  <IS0UT,5)  K,UKNL, VKNL,WFZ(K,N,LI , OXKNL,OY KNL , DFZ <K,N,L1 ,PXNL 

920  CONTINUE 
IF(MCC1I.NE.2»  GO  TC  9205 
WRITE (ISCUT, 61 

WRITE (ISOUT, 7)  <K, U$UM(K,N,L> ,VSUH<K,N,L> ,DXSUM( K,N?L1 , 

2DYSUM (K,N»L),RSUM(K,N,H»K=1,K8HXI 
9205  0<AV=0<AV/ZSPAN 

WRITE(ISOUT,8i  L,N,  DKAV 

CALL  CNTRtNET ,NETSU,N,XG,YG, ICF,JCF »NCF! 

XQ=XG 

Y  y  Q 

CALL  NEST<NET,NETSU,XQ,YQ,NOATG,XL,XR,YL,YU,ICF,JCF,NCFI 

AREAN=  (XR-XH  •  (YU-YLI 
OAVG (LI =DAVG(L I ♦OK A V* A RE AN 
00  9210  KL=1, KBHX 

9210  WAVG (<L, L) -  WAVG(<L,L)  ♦  WFZ ( «L,N, L) * AREAN 

921  CONTINUE 
OAVG(L)=OAVG(L}/AREA 
DO  9215  KL=1, K0HX 

9215  WAVG(KL,LI=WAVG<KL,L)  /  AREA 

WRITE (ISOUT ,151  L,  0AVG(U 

WRITE(IS0UT,16I 

WRITE ( ISOUT ,17) (<,W8VG(<,L> ,<=1,<9HX> 

922  CONTINUE 
RETURN 
END 


AOMTN117 

AOMIN118 

AOMIN119 

AOHIN120 

AOHIN121 

A0MIN122 

AOMIN123 

ADMIN124 

A0MIN125 

ADMIN126 

AOMTN127 

ADMTN128 

A0MIN129 

AOMIN130 

AOMIN131 

AOMIN132 

ADMIN133 

A0MIM134 

AOMTN135 

AOMIN136 

AOMTN137 

AOMIN138 

A0MTN139 

AOMIN140 

ADMIN141 

ADMTN14? 

ADMIN143 

A0MTN144 

A0MIN145 

A0MIN146 

AOMIN147 

A0MIN148 

ADMTN149 

ADWIN150 

ADMIN151 

AOMIN152 

A0MIN153 

AOMIN154 

AHM1N156 

ADMIN156 

A0MIN157 

AOMTN150 

APMIN159 

AOMIN160 

ADMIN161 

A0MIN162 

A0MIN163 

A0MIN164 


SUBROUTINE  AOVEC ( NE T,NETSU,ZBH, TI HUP, USUMvVSUM,DXSUM,DYSUM, RSUH, 
1T0EP,CAV,PMAS,PSIZ,ICF,JCF,NCF,KBHF,N0ATF,LTIMF,CAVS,WFZ) 

C  SEPTEMBER,  1971 

C  SUBROUTINE  AOVEC  TRANSPORTS  PARCELS  BY  SIPPLE  AOVECTION  °LUS 

C  SETTLING.  PARCEL  TOP  ANC  BASE  ARE  TRANSPORTED  SEPARATELY,  ANO  THE 

C  RESULTS  ARE  SHEARED.  THE  COMMON  VARIABLE  ZP  IS  REDEFINED  HEREIN. 

C  ZP  PARCEL  CENTER  Z  COORDINATE  BEFORE  AOVECTION,  ECXEPT  AS 

C  REOEFINED  IN  SUB.  AOVEC 


COMMON 

/QPARM' 

IPOUT 

, IPARIN, NBLK 

,NAT 

,NOELT 

,KX 

,  KKH 

l»NSEQO 

,  TCX 

,JCX 

, NCX  , K8HX 

,NDATX 

,LTIMX 

,  TSIN 

,  ISOUT 

2,ED0Y 

, FHBEL 

,LSTEP 

, MC(1S) , HINT 

,XLLC 

,  YLLC 

, THETA 

,ZH1’N 

IfCSKlP 

,MINT 

,  7MAX 

, TIMEX  , OT 

»DZ 

,XP 

,YP 

*ZP 

4,DINCR 

,  DOWN 

♦  TP 

, ZLOH  , DHAF 

,RHAF 

.ROPART 

,ZUPP 

,  VETA 

5.00PFN 

,  CROSS 

,TIME 

, KKMA1  , KKX 

, KKXS1 

, KXH IN 

,NOATP 

DIMENSION  NEK ICF, JCF) , NET SU  (NCF) , Z8H ( KBHF) ,USUH( KBHF, NCATF, LTIMF) 
DIMENSION  VSIJH  (KBHF  ,NDATF,LTIMF!t  ,  OXSUM  (KBHF, NO  ATF,  LTIMF) 

DIMENSION  DYSUM(<8hF,NDATF,LT IMF) , TIMUP (LTIMF) 

DIMENSION  RSUM(KBHF,NOATF, LTIMF! 

DIMENSION  CAVS(KBHF) ,HFZ( KBHF, NOATF, LTIMF) 

MC3=MC(3) 

EPS=0.1 

N0EP=0 

Z0EP=ZMIN 

CHANGE  ZP  FROM  PARCEL  CENTER  TO  PARCEL  BASE  ALTITUDE. 

7®=7L0W 

CALCULATE  TRANSPORT  OF  PARCEL  RASE. 

IF  (  (ZP-7DEP) «LE  .EPS*  GO  TO  1411 

CALL  TRANP (NET, NETSU,7BH, TIMUP, USUH,VSUM,DXSUM,DYSUH,oSUM, 

1NDEP,TOEP,ZOEP,XOL, YOL ,ZOL,TOL,SIGXL, SIGYL, ROL ,NDATL, ICF, JCF,NCF , 
2KBHF,N0ATF, LTIMF, 0,CAVS,WF7I 
GO  TO  1417 

1411  TOL=TP 
XOL=XP 
YOL=YP 
ZOL=ZP 
ROL=0. 

SIGXL=RHAF 

SIGYL=RWAF 

CHANGE  ZP  FROM  PARCEL  BASE  TO  PARCEL  TOP  ALTITUDE. 

1412  ZP=7L0W«-0WAF 

CALCULATE  TRANSPORT  OF  PARCEL  TOP. 

IF (  ZP-ZOEP  .LE.EPS)  GO  TO  1414 

CALL  TRANPfNET, NETSU,ZBH, TIMUP, USUH,VSUM, OXSUM, DYSUM,RSUM, 

1NOEP,TOEP,ZOEP,XOU,YOU,ZOU,TOU,SIGX»J,SIGYU,ROU,NOATU,ICF,  JCF,NCF, 

2KBHF, NO ATF, LTIMF, 0 , CAVS»HFZ) 

GO  TO  1415 

1414  TOU=TP 
XOU=XP 
YOU=YP 
ZOU=ZP 
ROU=0. 

SIGXU=RHAF 

SIGYU=RNAF 

CALCULATE  SHEAR  OF  PARCEL  TOP  AND  BASE  RESULTS. 

1415  ZOUTN- (ZOL^ZOUI  7  2. 

TOt«TN=(TOL*TOU)/2. 

IF  (ASS  (XOU-XOl.l .  GE.  1.0E-30)  GO  TO  1404 
IF < ABS (YOU- YOL) ,GE . i. OE-30)  GO  TO  1403 


ADVEC 

1  | 

ADVFC 

2  I 

AOVFC 

3  | 

ADVFC 

4  ! 

ADVEC 

5  | 

AOVFC 

6 

AOVEC 

7  i 

ADVEC 

8 

ADVFC 

9  ! 

ADVFC 

10  i 

ADVEC 

11  i 

ADVEC 

17 

AOVFC 

13 

Aovrr 

14 

ADVFC 

15 

AOVFC 

IF 

AOVFC 

17 

AOVFC 

13 

ADVEC 

19 

ADVEC 

20 

ADVEC 

21 

ADVFC 

22 

AOVFC 

23 

AOVFC 

24 

ADVFC 

25 

ADVEC 

26 

ADVFC 

27 

ADVFC 

26 

ADVEC 

29 

ADVEC 

30 

ADVEC 

31 

ADVEC 

32 

AOVEC 

33 

ADVFC 

34 

AOVEC 

35 

ADVEC 

36 

AOVEC 

37 

ADVEC 

33 

ADVEC 

39 

AOVFC 

40 

ADVEC 

41 

ADVEC 

42 

ADVEC 

43 

ADVEC 

44 

AOVEC 

45 

AOVEC 

46 

AOVEC 

47 

ADVEC 

48 

ADVEC 

49 

ADVEC 

50 

AOVEC 

51 

AOVFC 

52 

AOVEC 

53 

ADVEC 

54 

AOVEC 

55 

ADVEC 

56 

ADVEC 

57 

ADVEC 

58 

ROUTN=0. 

GO  TO  1*05 

1403  POUTN=l. 57070633 
GO  TO  1*05 

1404  ROUTN= AT  AN ( CYOU-YOL  )  / IXOU-XOLI 1  *  i 

1405  R=POUTN-ROL 

SXL  =  l./S3Rr(!C0S<RWSIGXL)**2«-CSINm  7SIGYL>*»2» 

SYL=1./SQRT  tCSIN(R)  /SIGXL)  **2MC0S  !R1/SIGYL>**2> 

R=POUTN-ROU 

SXU=i./SQRT<(COS(R> /SIGXU)  **2MSINtR>  7SIGYU>**2> 

SYU  =  i./SORT  CCSIN(R*,  7SIGXU)  **2MC05  (P»  /SIGYU>»*21 
SXOTN=  {SXU^SXL^QRT  {{XOU-XOLJ  **2M  YOU-YOLi  **2  >  172. 

SYOTN=SORT  «SY‘J*SYL» 

XOUTN=XOi.MSXOTU-SXtl*COS*RO>JfN'j 

YOUTN=YOLMSXOTN-5:XO*STN(ROUTN1 

CALL  CUMBER (XOL'TW* YCUfN» ZOUTNj TOUTN* SXOTN^SYOVN, P^ftS.PSIZ »ROUTN,  0 
1T50UT,  IPO’JT,HC3,NRLK) 

RETURN 

END 


AOVEC  50 
AOVEO.  60 
AOVEC  61 
AOVEC  62 
AOVEC  6? 
AOVEC  64 
AOVEC  65 
AOVEC  66 
AOVEC  67 
AOVEC  6P 
AOVEC  60 
ADVEC  70 
AOVEC'  .71 
AOVEC  7? 
AOVEC  73 
AOV.EC  74 
AOVEC  75 
AnyEr  76 
AOVEC  77 
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SUBROUTINE  AMBNT C ZCH,DF2,WFZ , ZHT,W ,0 IFF, A A, 99 , CC ,nEN0M,F , F, NET ,  AMPNT 

iN£TSU,ZPH,TIHUP,USUM,VSUM,OXSUM,DVSUM,RSUM,NTRIP,LTTM,ICF,JCF,NCF,  AMPNT 
ZKBHFjNCATF-LTINF, KKF,C*VSi  AMPNT 

SEPTEMBER  1971  AMPNT 

SU»0»JTINE  AMBNT  DETERMINES  FROM  INPUT  FALL  PATE  ANO  VERTICAL  KINO  AMPNT 
AND' TURBULENCE  DATA  THE  COEFFICIENT  DATA  VECTORS  AMPNT 

AA?!O=S2MOIFF(<+tt*0IFF(<))*Sl*F(<>l)  ,  AMENT 

BB(K)=S2»(0IEF  CKfl»  ♦2.*DIFF<t-  I  i-OIFF  (K- 1) »  f  S1*F  f  K)  ,  AMPNT 

CC(K1=S2*(0IFF|(O  +  0SFF(i<-15>  ,  AMPNT 

DEN0MU1=l*-THETA*{<3F'X)-CC(«O  *E  (K-11  )  ,  AMPNT 

E (X ) =THFTA* AA  tKJ/OFNOMtK) ,  AMENT 

WHERE  AMPNT 

•  Si=0T70Z,  AMPNT 

S2=0T/ (2.* iOZ) **2) ) ,  AMPNT 

ANO  F  IS  a  TEMPORARY  WORKING  SPACE  FOR  VERTICAL  VELOCITIES.  AMPNT 

NOTE  THAT  0(0  FOR  THE  NEXT  TIME  STEP  IS  GIVEN  BY  AMPNT 

A A  (K) *QCK*1>-9BIK) *CtK) *CC CK) *Q(K- II .  A-PNT 

DF<XS1-  VE°TTC AL  DffUSIVlTY  AT  ALT ITUOE  INC°EMENT  <X-l  AMPNT 

LT1M  -  A V M OS •  UPDATE  INDEX  FOR  ARRAYS  OFZ  ANO  Wf?  AMPNT 

NTRIP  -  OPTION  COOE  FOR  H0RI70NTAL  AOVECTION  OF  MASS  ALOFT  AMPNT 

POSITIVE  IF  NOATO  IS  STORED  IN  NTRIP  AMPNT 

NEGATIVE  IF  NOATO  IS  TO  BE  FOUND  Via  CALL  TO  SUB.  T°ANP  AmrnT 
COMMON  /QPARM/  1P0UY  , IPAR IN , N8L<  >NAT  sU0ELT  ,<X  ,<KM  AMPNT 

‘i,NSEOC  ,ICX  ,JC*  ,NCX  , KBHX  ,NDATX  ,LTIMX  ,TSIN  , ISOUT  AMENT 

' 2, EDDY  t FMBEL  fLSTtr  ,HCtl5),WINT  ,XLLC  , YL.LC  , THETA  ,ZMlN  AMPNT 

3,CCKIP  ,MINT  ,ZMAX  , TIMEX  ,0T  ,6?  ,XP  ,YP  , ZP  AKPNT 

4.0INCR  t DOWN  fTP  , ZLOW  ,OWAF  -PKAF  ,RDPART ,7UPP  >  VET  A  AMPNT 

5,D0PEN  , CROSS  ,  TIME  ,  k'KMAl  ,KKX  ,XKXS1  VKXMIN  ,NOATP  AMPNT 

COMMON  /QOPLE/OFKXS 1 ,£  FFLUX,  FMA8,PHI,THET0  AMPNT 

DOUBLE  PRECISION  A  A  ?XXF  I  ,98  (k~KF»  ,CC  (KKFI  ,CENOM  (XKF)  ,E  T  KKF)  ,F  CKKF  )  AMPNT 
DIMENSION  OIFF(KKF)  ,WtKKF)  ,7HT(«F)  ,OFZ  «<BHP,  NEAT  F,LT  IMF!  AMRNT 

DIMENSION  WFZ(KBHF,KOATF,LTIMF) ,TIMUP<LTIMF) , Z9H (KRHFT ,ZCH(KBHF)  ANBNT 
DIMENSION  axSUM(K9HF,NnATF,LTIHF» .OYSUMtKBHF, N3ATF,LTIMF)  AMRNT 

DIMENSION  USUMtKP>-F,NOATF,LTIMF}  ,  VSUM*K3HF,MDATF,LTIMF-,  AMPNT 

DIMENSION  RSUM(KBPF,NrATF,LTIMF) , NFTTICTtJCF» .NETSUtNCF)  AMPNT 

DIMENSION  CAVSiKBHF)  AMPNT 

DOUBLE  PRECISION  0FXXS1, EFFLUX, FMAB, PHI, THETQ  AMPNT 

DOUBLE  PRECISION  S1.S2  AMPNT 

DATA  PROGRM/E AHBNT  it  AMPNT 

CONSTRUCT  F  AND  OIFF  FOR  K=B,*..,KX  AMRNT 

NOEP=tOG  AMRNT 

TQHP=TIME  AMPNT 

NOATO=IABS<NTRIP|  AMRNT 

COMPUTE  KBH  ANC  INITIALIZE  AMPNT 

QIFFt<<M|=0.  AMPNT 

F(KKH)=C.  AMRNT 

ZOLO=ZHT(KKM*  AMPNT 

OOLD=OIFF<KKM»  AMPNT 

FOLO=FTKKM»  AMPNT 

CO  1  <=I ,<8HX  AMRNT 

KBH=K  AMBNT 

ZNEW=ZCHUBH)  AMBNT 

IFfZOLO.LT.ZNEWI  GO  TO  2  AMBNT 

1  IFtKBH.EQ. KBHXI  CALL  ERROR TPROGRM, -1 , ISOUT)  AMBNT 

2  lFtNTRIP.LT.-U  AMBNT 

1CALL  TRANPtNET,NETSU,ZBH,TIMUP,USUN,VSU«,DXSUH,OYSUH,RSUM,  AMPNT 

2N0EP»T0EP, ZNEW,70, YC,ZC,TO,SIGXO,SIGYO,RO»NOAT  0, ICF, JCF,NCF,KBHF,  AMPNT 
3W0ATF,LTIMF,l,CAVS,KF?t  AMRNT 
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18 
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5g*  £®S* 


nNEH=Of-Z(KRH;,NOA  TO.LflNI  AMBNT  69 

FNEW=-WF7<<RH,N0AT0  ,LTIM)  AMPNT  60 

OSLOPE=  inNFW-OOLD)  /  (ZNEW-70LD5  AMPNi*  63 

FSLGPE  =  (FNEW-FOLm  /  (ZNEW-70LDS  AMBNT  6 2 

DO  5  KX=KKMA1,KKX  AMPNI  63 

7HTKK=ZHT<KK»  AMPNT  64 

IFC7hT<<.lT.ZNEH)  GC  TO  4  AMPNT  66 

JF  (K8N.LT.KPKX-1I  GC  TO  3  AMPNT  66. 

40  DIFF(i<<|=ONEW  AMPNT  67 

F(KtO=FNEH  Amp  fit  6" 

GO  TO  6  AMPNT  6° 

3  KBH=K°M* 1  AMPNT  70 

IF(ZKT«.GE.7CH|<BHn  GO  TO  30  AMpNT  71 

QOLD=CNEW  AMOKT  7? 

FOLO=FNFW  AMPNT  73 

70LP=ZNEW  AMPNT  74 

GO  TO  38  AMPNT 

30  AMf-'-MT  76 

TFiKBC.L  T.XBKX-1)  GC  TO  ?,?  AMPNT  77 

<BH=<BC  AMPNT  78 

31  ZNEH=7CH(KBM1  AMPN?  79 

IPINTPTP.LT.-1  }  AMPNT  80 

1CALL  TRANP(NET»NFTSU»ZBH»TIMUP»USUMt VSUM  >DXSUM»QYSUM*RSUMf  AMPNT  81 

2NnFP,T0EP,ZNEU,XO, 7C,Z0«T0tSIGX0,SI6Y0eP0,NPAT0, ICF , JCF, NCF, <PHF ,  AMPNT  A? 
3NBATF, LTIMF,lf CAVS, 4FZ)  AMPNT  83 

ONEW  =  OFZ(K8H,NOATO,LTm  AHPN7  84 

FNfW--WF2t<RH,N0ATOfLTIM)  AMPNT  85 

GO  10  4G  AMBNT  86 

32  <RSXM1=KBH3-1  AMPNT  8 7 

DO  35  K=K3CtKBHXMl  AMPNT  38 

K9H  =  (C  AMPNT  89 

IFCZHTkK.GF.ZCHOCBH)!  GO  TO  35  AMPNT  90 

ZOLO=7CH{KRH-ll  AMPNT  91 

IFtNTPIF,tT.-l»  AMPNT  92 

1CALL  TRANPlNET,NETSU,ZeH,TrMUP,U?:L!M,VSiJM,bXSUM,nVSUMt»SUM,  AMPNT  93 

2NDEPrTOEP.ZOLO,XO,  YCf  ZC* TO ,5 IGMO, S  IGt 0 , POTl'!D£¥  Or.  IC5- ,  JCF,  NCPV  <RHF ,  flMpNT  94 
3NOATF,LTIMF,l,CAVS,i«FZ>  AMPNT  95 

COLG-rFZi<RH-l  fNDATC,tTTni-  AMP*}'  96 

F0LD^-WFZCY8H-1*N0ATC,LTIM)  AMPNT  97 

GO  TO  36  AMPNT  98 

35  CONTINUE  AMPNT  39 

CO  TC  31  AMPNilOO 

38  ZNEW=ZCH<<RHJ  AMBNT101 

I?(NTRIP.LT.-1*  AMeNT102 

1C ALL  TRANP{NET,NPTSU,ZBN,TrMUP,USUM,VSUM,OX^UMtOYSUM,RSUM,  AMPNT 10  3 

2«0EP,TDEP,ZNE«tX0, Y C * 2G. i O.SICXO, SIGYO,RO,NO«TO, ICF, JCF»NCFr KPHF ,  AMFNT1U4 
3NOATF,LTIKP,l,CAVS,VFZ»  AMBNT 105 

0NEW=0F7(K3H,N0AT0, LTIM1  AMPMT106 

FNEW=-HFZ(KBH;NDATO,l.TIHI  AMBNT1H7 

OSLOFE  = (ONEW-OOlH) /  (ZNEW-ZOLDl  AWPNT1GC 

FSLOPE  =  <FNEW-FOI.O)  x?ZNCW-ZOL01  AMBNT10Q 

4  0IFRi<O=DUt.0*-0Sl  0PF»17H"<*<-70L0)  AHBNTlltJ 

F  CKK? =PCLB6FSL0PEo T  ZHTXK-ZOLOI  AMSNTlli 

5  CONTINUE  AMPNT H2 

CORRECT  OIFFUSlVITItS  FOF  CROSSING  TRAJECTORIES  EFFECTS  AHPNT113 

00  150  «=<XM,<*X  ASPNT 114 

150  CIFFfKK»=00«N*i)T.PF(KX>  AHPNT115 

COMBINE  STILI.-AIR  FALL  RATES  WITH  VERTICAL  WINDS  AMPNT  116 


Wtf> 4: W ***  •&  . . .  ,»w« 


no  250  KK“KKM,KKX 

250 

COMMUTE  ECO)  IN.  ACCORDANCE  WITH  REFLECTIVITY  CONTROL  VARIABLE  VETA 
>  {XKM$=<J. 

IFCVETA.lt. fl.)  GO  TC  17 
HREF =  F  (KKM  Ail 

3RFF=CQIFFl«>mOIFF<«MAin/x2.*DZ> 

ECKKM1  -(WWEFfQRCF)  t  (YETA*9REF» 

COMPUTE  AA,B8,CC,0EN0K,  ANC  E  FOR  <=i , .  .  .  ,KX-i 
17  S2  =  07 
S1  =  0T 
S1=S1/S2 
C2=S1/(2.*S2> 

00  18  «=<<MAl,KKX 
IS  CCCKf<)=?2*  COIFFtKKl  *OTFF<KK-in 
00  IB  «=K<M,<,<XS1 

19  AA(KtO=CC(KK*l>fSl*F(KK  +  li 
00  20  «=<<MH,KKXS1 

96 ( KK) =  AA(KK-1)«-CC (XKf II 

DEN0M('lO=i.»-THETQ*  (88 CKO  -CCCKO  *E(KK-1I ) 

20  E(KKl=THETG*AA(KK>  7CENCMCKK) 

0FKXS‘1~DIFF  («XS1) 

C  HRITE(IS0UT,2221» 

C2221  FORMAT  (<OCONTENTS  OF  ARRAYS  ZHT,  F,  DIFF,  AA,  BP,  CC,  OENOM,  E*> 
C  WRITE (ISOUT .2222)  Z  FT 

C  WRITE ( IS OUT ,2222)  r 

C  WRITE CISOUT, 2222)  DIFF 

C  WRIT  F i ISOUT ,2222)  AA 

C  WRITE  f  ISOljTt2222)  9P 

C  WRITE <ISOUTf222t)  CC 

C  WRITEf ISOUT, 2222)  OENOM 

C  HRITECISOUTs2222:  E 

C2222  FORMAT ^ #0* < 13E 10. 2 ) ) 

RETURN 

END 
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SUBROUTINE  SDRVS(NET,NETSU,ZBH,ZCH , TIMUP, USUM,  VSUM,OXSUM»  OYSUH  ,  SPPVS 
1RSUM,0FZ,WF2,GAVG,WAVG,ALT,RH0,ATEMP,AA,8B,CC,DEN0M,0IFF,£,F,G,W,  SpPVS 
2ZH  T, ICF, JCF,NCF,KPHF,NCATF,LTIMF,KKF,NATF,CAVS)  SPPVS 

SEPTEMBER  1971  SPPVS 

SUBROUTINE  SPRVS  SUPERVISES  DIFFUSIVE  AND/OR  ADVECTIVE  TRANSPORT  SPPVS 
OF  FALLOUT  PARCELS  LISTFO  ON  TAPE  IPARIN.  PARCEL  PARAMETERS  ARE  SPPVS 
STORED  IN  ARRAYS  XPAR, YPAR,ZPAR, TPAR, P0AM,PSAM ,RWFR,DWFR, ZLWF, VWFRSPRVS 
ONLY  ONE  PARCEL  IS  TRANSPORTED  AT  A  TIME.  FOR  THIS  PARCEL  ABOVE  SPPVS 
ITEMS  ARE  STORED  IN  XP, YP, ZP, TP, PS IZ.PMAS ,RKAF ,OWAF, ZLOW, VWAF.  SPPVS 

XPAR  -  X  COORDINATE  OF  PARCEL  CENTER  DATA  VECTOR  CAT  TIME  TPAR)  SPPVS 

YP AP  -  Y  COORDINATE  OF  PARCEL  CENTFR  DATA  VECTOR  (AT  TIME  TPAR)  SPPVS 

ZP  AR  -  Z  COORDINATE  OF  PARCEL  CENTER  DATA  VECTOR  (AT  TIME  TPAR)  SPRVS 

TPAR  -  TIME  OF  DEFINITION  OF  CLOUD  PAPCEL  DATA  VFCTOR  SPPVS 

PD AM  -  MIDPOINT  OF  PARCEL  PARTICLE  SIZE  CLASS  DATA  VECTOR  SPRVS 

PS AM  -  TOTAL  MASS  CF  PARCEL  DATA  VlCTQP  (AT  TIME  TPAR)  SPPVS 

OWPR  -  RADIUS  OF  PARCEL  AT  C.  0.  M.  DATA  VFCTOR  (AT  TIMF  TPAR)  SDPVS 

OWFR  -  PARCEL  THICKNESS  OATA  VECTOR  (AT  TIME  TPAR)  SPPVS 

ZLWF  -  ALTITUDE  OF  PARCEL  8ASA  DATA  VECTOR  (AT  TTNE  TPAR)  SpPVS 

VWFR  -  PARCEL  VOLUME  DATA  VECTOR  (AT  TIME  TPAR)  SPPVS 


COMMON 

/QPARM / 

IPOUT 

♦IPARIN 

,  NBLK 

,NAT 

,NDELT 

,<X 

,KKM 

SORVS 

l,NSEQO 

,rcx 

,  JCX 

,  NCX 

,  KBHX 

,NOATX 

,LTIMX 

» ISIN 

y ISOUT 

SPPVS 

2*E09Y 

, FMBEL 

,LSTEF 

,  MC (1 3) 

,  W INT 

,  XLLC 

»  YLLC 

, theta 

,  ZMIN 

SpF  VS 

3,CSKIP 

,  MINT 

,ZMAX 

,  TIMEX 

» OT 

,DZ 

,XP 

,YP 

,ZP 

S°PVS 

4, OTNCR 

♦  DOWN 

» TP 

y  ZLOW 

,  OWAF 

,RWAp 

, RO° ART 

,ZUPP 

,  VFT  A 

SPPVS 

5 , OOPEN 

♦CROSS 

» TI  E 

,KKMA1 

,«X 

, KKXSl 

,<XMIN 

,NOATP 

SPRVS 

COMMON  /Q09LE/0FKXS1, EFFLUX, FMAB, PHI, THETC  SPPVS 

DOUBLE  PRECISION  A  A  (KKF)  ,88 (KKF) , CC  (KKF)  ,OENOS  <<KF)  ,E  («F> , F  (<<F )  S°PVS 
DOUBLE  PRECISION  0 ( KKF)  SPPVS 

DIMENSION  ALT(NATF)  ,RHO(NATF) ,ATEMP(NATF)  SPPVS 

DIMENSION  NET ( ICF, JCF) , NETSU(NCF) , ZBH (KBHF) ,ZCH( K8HF) ,DIFF (KKF)  SPRVS 
DIMENSION  USUM ( KBHF ,N0  ATF  ,LT  IMF)  ,  V  SUM  (KBHF,  NDA  TF ,LTIMF) , W (KKF )  SPRVS 

DIMENSION  DXSUM(<BHF,N0ATF,LTIHF),DYSUM(<9HF,NDATF,LTIMF>  SpPVS 

DIMENSION  0FZ(K8HF,N0ATF,LTIMF) ,ZHT(KKF)  SPPVS 

DIMENSION  WFZ ( KBHF , KOATF,LT IMF) ,DAVGILTIMF) ,WAVG (KBHF, LT IMF)  SPRVS 

DIMENSION  TIMUPCLTIMF) ,RSUM ( KBHF, NDATF, LTIMF)  SPPVS 

DIMENSION  XPAR(IOO)  ,YPAP(100) ,ZPAR (1001, TP AP{ 100) ,PDAM(100>  SPPVS 

DIMENSION  PSAMC100)  ,RWFR (100 ) ,OWFR f 1 0 0 ) , ZLWF( 1 00 ) ,VWFR(100)  SPRVS 

DIMENSION  CAVS(KBHF)  SpPVS 

DOUBLE  PRECISION  OFKXS 1 , EFFLUX, FMA 8 , PHI , THETQ  SPRVS 

8014  FORMATf#«-#T102,E12.4,I4)  SPRVS 

8015  FORMAT(#+#T103, /AIRBORNE  (ADVCN)/)  SPRVS 

3016  FOPMAT (/  #14, 3E12.4)  SPPVS 

8017  FOPMAT(/*/T103»/AIRf?0RNF  ( OIFFN) /)  SPRVS 

3018  FORMAT (/4-/T103»/ADVECTIVE  TRNSPT/)  SPRVS 

8019  FORMAT(#«-/T103,/  IMPACTED  WAFER#!  SPPVS 

3020  FORMATf/+/T103, /OUTSIDE  WINDGPIO/)  SPPVS 

8021  FORMAT(#0#36X, /PARTICLE  SIZE  CLASS/E12. 4, /  MICRONS#)  SPRVS 

3022  FORMAT (#  #22X,#FALL  RATF#E12.4,#  METFRS/SEC  AT  ALTITUOE#E12. 4, #  SpPVS 

1  HETERS//22X, /UPPER  LIMIT  INPUT  ALTITUDE  FOR  ADV.  TRANSP.  IS/E 12 . 4SPRVS 
2,#  METERS#)  SPPVS 

8024  FORMAT (#0#T2, #NSEQ# Til, #XP#T23,#YP#T35,#ZP#T47,#TP#T5 8, /PMAS/T70 , /SPRVS 

1RWAF#T32,#ZL0W#T94, #OWAF#T107,#nZ#T115,#KX#/)  SPPVS 

8025  FORMAT (/ONEGAT IVE  OEPOSIT.  WAFER  NO./I4,#  AT  TIME/E12.4, #.  VARIABLSPRVS 

5ES  EFFLUX, FMAB,OEP  #2D12.4,E12.4/)  SPPVS 

DATA  PROGRM//SPRVS  #/  SPPVS 


JF=100 
K<MA1  =  «M«-1 
THFTQ=THET  A 
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PHI=1.-THETG 

KX8F=0 

MC3  =  MC  f 3 1 

nseo=o 

PSZBE=-2.0 

COMPUTE  OVERALL  AVERAGE  VERTICAL  VELOCITY  FOR  THE  FIRST  UPDATE 
WAVGK=Q.O 
KBHMl=KBHX-i 
DO  50  K=l, KBHH1 

50  WAVGK=WAVGK  f  WAVG « K, 11  * <ZBH (K*ll  -  ZBHfOl 
WAVGK=W«VGK/<ZBH<KOHX)-ZeH(l)  I 

COMPUTE  TIMEX  MARGIN  FACTOR  FOP  AOVECTIVE  TRANSPORT  AIRBORNE  TEST 
IF(NDATX-1)70,70,60 
60  SLOP=l.l 
GO  TO  30 
70  SLOOri.0 

COMPUTE  MINIMUM  SMALL  ALT ITUDF  INCREMENT  OZMIN 
80  OZMIN=(ZMAX-ZMIN»/KX 

Ci'E  I °ARIN  TAPE  AT  BEGINNING  OF  INPUT  PARCEL  RLOCK 
100  READ ( IPARINl  NP 

IFTNP.LE.0I  GO  TO  806 

IF(NP.GT.JF1  CALL  EPROR<PROGPH,-1O0,ISOUT1 
COPY  IN  A  8L0C<  OF  INPUT  PARCEL  PARAMETERS  FROM  TAPE  IPARIN 

REAOCIPAPTNI  f  XPAR  {  J) ,  YPAR  ( J)  ,  ZPAR  f  J1  ,  TPAP1  Jtl  ,  PD  AM?  j| ,»SAHf J1 , 
lRWFR(J) ,OWFR!Jl ,ZLWF(J1 ,VWFR? J1 ,J=1,NP) 

COMMENCE  PROCESSING  BLOCK  OF  INPUT  PARCELS  ONF  AT  A  TIME 
PO  1000  J=1,NP 
4SEQ=NSEQH 

IF (NSE0.LT •  NSEQOl  GC  TO  1000 
XP=XPAR«J) 

YP=YPA9?J1 
ZP=ZPAR(J1 
TP=TPARf J» 

PSIZ=1.0E6*POAMTJ) 

PMAS=PSAM(JI 
RWAF=RWFRf J>/2. 

DWAF=OWFR?Jl 

ZLOW=ZLWF?J> 

VWAF=VWFRf J1 

CHECK  FOR  NEW  PARTICLE  SIZE  CLASS 

IFTABSUPSIZ-PSZBEI  /PSI7) , LE. 1 . 0E- 101  GO  TO  103 
WRITE  CISOUT ,8021)  PSI? 

COMPUTE  MID-ATMOSPHERE  FALL  RATE  FAV  FOR  NEW  PARTICLE  SI7E  CLASS 
H=CZMIN+ZMAX1/2. 

CALL  TRPL(H, NAT, ALT, RHC, DENI 
CALL  TRPLf H, NAT, ALT, ATEMP, VIS! 

CALL  FALRT(P5IZ,R0FART,H,0EN,VIS,FAV,IS0UT> 

FAV-FA V- WAV GK 

COMPUTE  UPPER  LIMIT  ALTITUDE  FOR  AOVECTIVE  TRANSPORT  OF  THIS  SIZE  PART. 
CALL  CALIBCZBH,KBHX,ZMIN,-1,KBHZI 
CALL  TRPLCZBHCKBHZ) ,NAT, ALT, ATEMP, VIST 
CALL  TRPLCZBH ( KBHZ1 ,NAT,ALT, RHO,DEN) 

CALLFALRT (PSIZ ,ROPART *ZBHIKBHZl*DEN,VIS»CAV»ISOUTl 

TDEP=TP  *  <ZBHf<8HZ)  -  ZMINI  /(CAV  -  WAVG CKBHZ-1 ,1) I 

KBHH1=KBHX-1 

OO  1001  IZ=KBHZ,KBHM1 

CALL  TRPLfJCHC  IZ  I  ,NAT, ALT, ATEMP, VISI 

CALL  TRPL (ZCH(  IZ  I  , NAT, ALT, RHO, DENI 
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CALLEALRT ( PSIZ , ROP ART ,  2CH(  IZ  ) ,OEN, VIS, CAV, ISOUT) 

TOEP=TOEP  f  CZBHdZd)  -  78H f TZ) ) / (CAV-  WAVG(I7,i)1 
IF(TOEP.GT;.SLOP*TIMEX)  GO  TO  1002 

1001  CONTINUE 
ZLtM=5.0£4 
GO  TO  1003 

1002  ZLTM  =  78H(IZ«-1I 

1003  WPITE  (ISOUT, 8022)  FCV,H,7LIM 
WRITE(ISOUT,3024> 

COMPUTE  PARTICLE  FALL  PATE  TABLE  FOP  EACH  ATMOSPHERIC  STRATUM  WHEN 
C  A  NFW  PARTICLE  SIZE  IS  ENCOUNTERED 
DO  101  <.<Z  =  1,<BHX 

CALL  T  RPL ( ZCH { KKZ ) , NAI, ALT  , A  TEMP, VIS) 

CALL  T PPL (ZCH (<<Z) , NAT , ALT ,RHO,DENI 

101  CALL  FALRT (PS  17 , ROP AR T , ZCH (KKZ) ,OE N, V I S , e AVS ( KKZ ) ,ISOUD 
COMPUTE  OTFFUSIVITY  CORPECTICNS  FOP  NEW  PARTICLE  SIZE  CLASS 
n0WN=(FAV*EDDYI **2 
CPOSP=l.  /"SORT  (1.4  4.  *OOWN» 

DOHN=l  ./S1RT<  1  .»-OOWN> 

PSZCE=PSIZ 

103  WRITE(ISOUT,8016)  NSEO , XP, YP , ZP, TP , PMA S , RWAF, ZLOW,DKA F 
CANCEL  PROCESSING  OF  PARCEL  IF  IT  HAS  ALREADY  IMPACTED 
IF  (IFIX(OWAF).GT.O)  GO  TO  1200 
WRITEdSOUT, 80191 

CALL  nUMPER(XFfYP,7P,TP,  RWAF,  PWAF,  PMAS, PSIZ, 0 . , 0 , 

1ISOUT, IP0UT,MC3,NBLK» 

GO  TO  1000 

COMPUTE  INOFX  OF  MESH  OR  SUB-MESH  CONTAINING  PARCEL  CENTER  POSITION 

1200  CALL  NEST(NET,NETSU,XP,YP,NOATP,XL,XR,tL,YU,ICF, JCF,NCF» 
CANCEL  PROCESSING  OF  PARCFL  IF  IT  IS  INPUT  OUTSIDE  ATMOS. 
IF(NOATP.GT.0»  GO  TC  1248 
WRITE ( ISOUT ,30201 
GO  TO  1000 

COMPUTE  AVERAGE  FALL  RATE  CAV 

1248  ZUPP=ZLOW*DWAF 
ZL0=ZLCW-ZMIN 
ZUP=ZUPP-ZMIN 
H=ZMlN4-ZUP/2. 

CALL  TRPL(H,NAT,ALT  ,RHO,OEN) 

CALL  TRPL(H,NAT,ALT  ,ATEMP,VIS) 

CALL  FALRT(PSIZ,ROPART,H,OENs, VIS, CAV,  ISOUT) 

CANCEL  PROCESSING  OF  PARCEL  IF  IT  WILL  REMAIN  AIRBORNE  «Y  DIFFUSION 
OAV=OOWN*OAVG(i> 

CA V=CA V- WA  VG< 

TFLY=TIMEX-TP 

IFdFLY.LE.Q. )  GO  TC  1249 

CALL  ESTHCZUP,CAV,OAV,TFLV,PUP) 

CALL  ESTM(ZLO,CAV,OAV,TFLY,PLO) 

POEST=i.C- (PUP-PLO) /(ZUP-ZLO) 

IF (CSKIP.LE.POEST)  GO  TO  1250 

1249  WRITEdSOUT, 3017) 

GO  TO  1000 

COMPUTE  TRANSPORT  BY  ADVECTION  IF  TRUNCATION  ERROR  IS  EXCEEOED 

1250  IF(OZMIN.LE.0.2*DAV/FAVI  GO  TO  1500 

CANCEL  PROCESSING  OF  PARCEL  IF  IT  WILL  REMAIN  AIRBORNE  BY  ADVECTION 
IF(ZLOW.LT.ZLIH)  GO  TO  1409 
WRITE (ISOUT, 3015) 

GO  .0  1000 
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SPPVS157 
SPRVS158 
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SPRVS166 
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1409  HRITEtTSOUT, 80131 

CALL  AOVECCNE  t*N£TSU,Z8H,TIMUP,USUM,VSUH,DXSUM9DYSUH»RSUM, 

1T0EP,CAV,PMAS,PSI7, ICF, JCF,NCF, <3HF, NDATF,LTIH~, CAVS, WFZ» 

GO  TO  1000 

COMPUTE  SMALL  ALTITUDE  INCREMENT  OZ 
1500  OZMAX=OWAF 

GO  TO  1502 


1501 


IF(OZMAX.GE.G2MINJ 
OZMAX=2.0*DZNAX 
GO  TO  1501 


1502  f)7=0*2*OAV/FAV 

IF(DZ,GT.3ZHAXI  D7  =  rZM.*X 
KX=iZ*AX-ZMINI/UZ4-l. 

IFCXX.LT.KXMTNT  KX  =  KXMlN 
07=IZMAX-ZMIN*  /KX 
WPTTEdSOUT  ,30 IV)  02, XX 

COMPUTE  ALTITUDE  INC  RE HF NT  AND  FALL  RATE  DATA  VECTORS 
IF(<-  .EO.KXPEl  GO  TC  1305 
KKX-^KM+KX 
KKXS1=KKX-1 
00  130V  KK=l,Kk'X 
ZHT<KK1=ZMIN+D7* (KK-KKM1 
CALL  TRPLtZHTtKO  ,  N  AT  ,  ALT,  RHO  ,OEN» 

CALL  TRPL(ZHTCKK» ,NAT, AL1 , AVEMP,VIS» 

130V  CALL  FALRT (PSI7,R0FAPT,ZHTI<<I ,DEN,VIS,M(KK)  , IS OUT? 

KXBE=KX 

COMPOSE  INITIAL  CONCENTRATION  DATA  VECTOR 
1305  CALL  CONCT7HT  ,Q,KKFT 

CANCEL  PROCESSING  OF  WAFER  IF  INITIAL  AIRPOPNE  MASS  IS  INADEQUATE 
IF(FMAg.LT.l)  GO  TO  1000 
DEP8E=0. 

NTRI?=NDATP 

COMPUTE  FIRST  DEPOSIT  INCREMENT  TIME  INTERVAL  TOELT 
TOELT=TFLY/NDELT 
TLARG=  TOELT 
TAE=V.*ZUP/!3. *CAV» 

IF (TAE.LT. TFLYT  TDELT= TAE/NOELT 
TSHAL-MINT*DT 

IF(TOELT.LT.TSMAL)  tdelt=tsmal 

COMMENCE  DEPOSIT  TIME  LOOP 
TIME=TP 
TPAUS=TIME 
LTIM=-1 
NPASS=1 
XOBE=XP 
YOPE=YP 

1  LSTEP=TOELT/DT  H. 

TOELT=LSTEP*OT 
135  TPAUS=TPAUS»TOELT 

COMPUTE  DATA  SET  TIME  INDEX  LTIM 

CALL  CALI9«tIMUP,LT!MX,TlHE,*l,LTIMA) 

IFaTIM.NE.LTIMAl  GO  TO  3 
IF:iA8S<NTPIP)-l»  31,32,31 
3  LTIM=LTIHA 

COMPOSE  CONCENTRATION  COEFFICIENT  DATA  ARRAYS 

31  CALL  AM8NT  f ZCH»DFZ,WFZ ,ZHT,W,0IFF,AA,B6»CC» OENGM,E »F,NET , 

1NET$U»ZBH»  TIKUP»USUM»  VSUM, OX SUM, OY SUM, R SUM, NT RIP , LTIM , ICF ,JCF, NC 
2KBHF,N0ATF,LTIMF,KKF,CAVS* 

NTRIP=-NOATX 


ts 
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SPPVS206 
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SPRVS209 
SPRVS210 
SPRVS211 
SPRVS212 
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SPPVS21V 
SPRVS215 
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SPRVS21R 
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SPRVS227 
SPRVS228 
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F , SPPVS230 
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COMPUTE  NUMfiEP  OF  ITERATIONS  ESTEP  OF  SOLUTION  TO  VERTICAL  DIFFUSION  SPPVS23? 
C  DIFFERENCE  EQUATION  FOR  THIS  DEPOSIT  INCREMENT  TIME  INTERVAL  SPPVS234 

32  IC(LTIM»GE.LTIMX)  GC  TO  138  SPRVS23F 

IFtTPAUS.LT «TIMUP(LTIM+1) )  GO  TO  140  SPRVS236 

TPAUS=TIMUP(LTIM*1J  SPRVS237 

TDELT=TPAUS-TIHE  SPPVS238 

LSTEP=TOELT/Of »i.  SPPVS239 

TOELT=LSTEP*DT  SPPVS240 

138  IFTTPAUS.LT. TIMEXI  GO  TO  140  SPPVS241 

TOAUS=TIHEX  SPPVS242 

TOELT=TPAUS-TIME  SPPVS243 

LSTEP  =  TDELT/'OT+l.  SPPVS244 

TOFLT=LSTEP*OT  SppVS?45 

IF(LSTEP.LE.O)  GO  TC  1000  SPPVS246 

140  IF (NPASS.NE.l)  GO  TC  4  SPPVS247 

NPASS=2  SPRVS248 

TQELTP=IOELT  SPRVS24? 

COMPUTE  DEPOSIT  INCREMENT  FRACTIONAL  MASS  DEP  SP*>VS250 

4  CALL  GIFFFCO,  AA,BP,CC,DENOM,E,FtOEP,«F)  SPPVS251 

CHECK  DEPOSIT  INCREMENT  FRACTIONAL  MASS  DEP  AGAINST  OOPEN  SPPVS252 

IF(DEP.GE.-DOPEN*TOELT)  GO  TO  5  SPRVS253 

HRITE(ISOUT,S025I  NSEO, TIME, EFFLUX ,FMAB ,OEp  SPpVS254 

5  IP (DEP.GE«DODFN*TDELT)  GO  TO  7  SDpVS255 

CHECK  CUMULATIVE  AIRBORNE  FRACTIONAL  MASS  FMA8  AGAINST  FM8EL  SPRVS256 

IF (SNGL (FMA3) . GT. FMPEL /  GO  TO  135  SPPVS257 

GO  TO  1000  SPPVS258 

COMPARE  PATE  OF  CHANGE  OF  DEPOSITION  RATE  DPCP  WITH  O.INCR  AND  ADJUST  SPPVS259 
C  NEXT  TOELT  SPRVS260 

7  OPOP=<COEP/ TUEin-CCEPBE/TOELTBII/rCELT  SPRVS261 

IF (OPOP. LT.OINCR)  GC  TO  IQ  SPRVS262 

TOELTB=TDELT  SPRVS263 

TOELT=TOELT/2.  S°RVS264 

IP(TOELT.LT.TSHAL)  TOELT=TSMflL  SPPVS265 

GO  TO  13  SPRVS266 

10  IFCDPDP.GT.OINCRt  GC  TO  13  S»PVS267 

TDELTR=TDELT  SPRVS268 

TOELT=2.*TnELT  SPPVS269 

IF(TOELT.GT.TLARG)  70ELT=TLARG  SPPVS270 

COMPUTE  DEPOSIT  INCREMENT  MASS  PHOEP  SPRVS271 

13  PMDEP=PMCS*OEP  SPRVS272 

OEPBE=L-EP  SPPVS273 

NDEP=0  SPPVS274 

ZOEP=ZHIN  SPRVS275 

TDEP=TIME  SPRVS276 

COMPUTE  DEPOSIT  INCREMENT  POSITION  (XO,YO,ZO»  AND  HORIZONTAL  DISPERSION  SPPVS277 
C  PARAMETERS  TStGXC , S IGYC ,R01  AT  TTME=TO  SPPVS27 8 

CALL  TRANP(NET,NETSU,Z9H,TIHUP,USUM,VSUM,DXSUMtOrSUM,PSUM,  SFRVS279 

lNnEPsTDEP»?OEPfXO,YC,ZO,TO,SIGXO,SIGrO,RC>NDATO,ICF,JCF,NCF,KBHF,  L‘PRVS20O 
2NDATF,LTIMF,1 , CAVS  ,  WF7J  SPRVS281 

CONTINUE  ON  TO  NEXT  WAFER  IF  THIS  ONE  LIES  OUTSIOE  WlNDFIELD  5°RVS282 

IF (NDATO.LE.O)  GO  TC  1000  SPRVS283 

COLLECT  FINAL  RESULTS  FOR  THIS  WAFER  AND  STORE  IN  BUFFER  DATA  VECTORS  SPRVS254 
XH=(XO*XOBE»/2*  SPRVS285 

YH= (YO*YUBET/2.  SP°VS296 

CALL  DUMPERTXHfYH,ZO,TO,SIGXO,SlGYO,PMDEP,PSIZ,RO,0,  SPPVS287 

1ISOUT, rP0UT,HC3,NBLK)  SPPVS288 

XOBE=XO  SPRVS289 

YOBE-YO  SPRVS29C 
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IF<TIME»DT,LToTIHEXI  go  TO  1 
1003  CONTINUE 
GO  TO  100 

COPY  OUT  BUFFER  OATA  VECTORS.  WAFER  PROCESSING  HAS  BEEN  COMPLETED 
806  CALL  QUHPERtO-,,0.,0.,0.,  0.,  0.,  0.,  0.,0.,999, 

1ISOUT,  IPOUT  ,NC3,N8LiO 

CALL  DUMPERfO.»0.»0.»0«,  0.,  0.,  0.,  Q.,C.,999, 

1ICOUT, IPOUT ,MC3,N8L<» 

REWIND  IPARIN 
ENO  FILE  IPOUT 
REWINO  IPOUT 
RETURN 
ENO 
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SUBROUTINE  TRANP(NET.NETSU,Z«iHf  TlrtUP, ’JSUH,VSUH,DXSUH,OYSUH,RSUM,  TRAMP  1 
1N0EP, TOEP»ZO£P,XO,YC,ZC, TO, SIGXO, SIGY0,R0,NDAT0,TCF, JCF, NCF, KBHF ,  TRANP  2 
2NOATF,LTIMF,<RIP,CAVS»WFZI  TRANP  3 

SEPTEMBER  1371  TPANP  4 

SUBROUTINE  TRANP  DETERMINES  « AT  AN  INPUT  TERMINAL  ALTITUDE  PLANE)  TRANP  5 
THE  WAFER  HORIZONTAL  CENTFO  POSITION  AND  DISPERSION  PARAMETERS  F  OR  TRANP  6 
AN  INPUT  TRANSPORT  FLIGHT  TINE  TRANP  7 

CAVS  -  PARTICLE  FALL  RATE  TABULATED  FOR  EACH  ALTITUDE  STRATUM  TRAMP  0 
TURBULENCE  X  COMPONENT  <WE IGHTEO  SUM)  3-OTM.  DATA  ARRAY  TRANP  9 

TURBULENCE  Y  COMPONENT  {WEIGHTED  SUM)  3-DIM.  OATA  ARRAY  TRANP  1C 

HORIZONTAL  SPACE  CONTROL  NET  MESH  2-DIN.  ARRAY  TRAN°  11 

HORIZONTAL  SPACE  CONTROL  NET  SUB-MESH  DATA  VECTOR  TD ANP  12 

WIND  HEAOING  ORIENTATION  ANGLE  (WEIGNTED  SUM)  3-OlM.  ARRAYTPANP  13 
WINO  X  COMPONENT  (WEIGHTED  SUM)  3-DIM.  DATA  ARRAY  TRANP  14 

WIND  Y  COMPONENT  (WEIGHTED  SUM)  3-OIM.  DATA  ARRAY  TPANP  15 

ATM0SoHERE  STRATA  BASE-ALT ITUDE  DATA  VECTOR  TD ANP  If 

COMPUTATION  MODE  SWITCH  TRAMP  17 

0  RAPIO  COMPUTATION  TPANP  18 

1  LAYERWISE  COMPUTATION  )®ANP  19 

OPTION  CONTROL  VARIABLE  TRAMP  20 

ZERO  IF  SIGXO  ANO  SIGYO  ARE  TO  BE  COMPUTED  TPANP  21 

NON-ZERO  IF  SIGXO  AND  SIGYO  ARE  NOT  TO  BE  COMPUTED  AND  TPANP  22 
IF  NDATO  IS  TO  BE  POSITIVE  ALWAYS  TRANP  23 

TOEP  -  AOVECTIVE  TRANSPORT  TIME  INTERVAL  TRANP  24 

KRIP  -  CONTROL  VARIABLE  TRANP  25 

0  FOR  AOVECTIVE  TRANSPORT  TPANP  26 

1  FOR  DIFFUSIVE  TRANSPORT  TPANP  27 

WFZ  -  VERTICAL  WIND  FIELD  TRANP  23 

ZOEP  -  AOVECTIVE  TRANSPORT  TERMINAL  ALTITUDE  TRANP  29 

TO  -  TIME  AFTER  PARCEL  ADVECTION  TPANP  30 

XO  -  PARCEL  CENTER  X  COORDINATE  AFTER  ADVECTION  TffflNP  31 

YO  -  PARCEL  CENTER  Y  COORDINATE  AFTER  ADVECTION  TRANP  32 

ZO  -  PARCEL  CENTER  Z  COORDINATE  AFTER  ADVECTION  TRANP  33 

SIGXO  -  PARCEL  MASS  HOR.  STANO.  DEV.  DOWNWIND  AFTER  ADVECTION  TRANP  34 

SIGYO  -  PARCEL  MASS  HOR.  ST4N0.  OEV.  CROSSWIND  AFTER  ADVECTION  TRANP  35 

NOATO  -  HORIZONTAL  iNOEX  OF  LATTICE  CELL  CONTAINING  POINT  (XG,YC>  TRANP  36 

RO  -  WIND  HEADING  ORIENTATION  ANGLE  AFTER  AOVFCTION  TRANP  37 


OX'SUM 

DYSUM 

NET 

NETSU 

RSUH 

USUM 

VS'JM 

Z3H 

MODE 


NDEP 


DIMENSION  NErtlCF,JCF)  ,NFTSU(NCF)  *  ZCH'KBHFV  ,  USUM  iK8HF,NDATF,LTlHFI  TPANP  44 
OIMENSIDN  VSUM(f'BHF,NDATF,LTIMP),OXSUM(<BHF,NOATE,LTXM?)  TRAN?  45 

DIMENSION  OYSUKCKBHF,NnATr  ,LTIMF*  „  TIMUR  (LUKE)  TRANP  46 

DIMENSION  CAVS(KBHF)  TRAP*  47 

DIMENSION  RSUM«<BHF,NOATF,LTIHF|  TRANP  40 

DIMENSION  WFZ(KBHF,NQATF,LTIMF)  TRANP  49 

DATA  PROGRM/*TRANP  7/  TRANP  50 

2  FORMAT  (7  TTME=7E12.4,7.  ALT=7E12. 4 , X-P0S=7E12.4t  7.  Y-P0S=7E12 .4TRANP  51 

2.7.  CELL=7I5,7  REACHED7)  TRAM®  52 

3  FORMAT (7  TIME=7Ei2. 4,7.  ALT=7E12.4„7.  X-P0S=7E12.4, 7.  Y-P0S=7E 12 .4TRANP  53 

3.7.  CELL=7I6,7  ATTEMPTED*)  TRANP  54 

4  FORD  AT  170  WAFER  WITH  INITIAL  CONFIGURATION  XP,YP,ZP,TP  74E12.4/r7REQTRANP  55 

10UIRE0  CHANNELLING  AT  CONFIGURATION  XC,YC,ZC,TC  74E12.4)  TRANP  56 

ERSILO=.0005  TRANP  57 

EPS=EPSlLO*WINT  TRANP  58 


COMMON 

/OPARM/ 

IFOU7 

, IPAp IN , NBL< 

,NAT 

-NDELT 

,<X 

,<KM 

TRANP 

38 

i,RSE0C 

,ICX 

,  JCX 

,NCX  , KBHX 

, NO ATX 

,LTIMX 

,  ISIN 

,  ISOUT 

TRANP 

39 

s 

l 

\ 

2 , EOOY 

5  FMCEL 

9LSTEP 

, MC (13) , WINT 

,xll:; 

,  YLLC 

, THETA 

,  ?MIN 

TRANP 

40 

3,CS<JP 

f  HINT 

,2MOX 

, TIMEX  ,DT 

,oz 

,x° 

,YP 

,ZP 

TRANP 

41 

j 

\ 

4,CINCP 

,DCWN 

.TP 

, ZLOK  , OKAF 

»RWAF 

, ROPART 

,Z  UPP 

,  VET  A 

TRANP 

42 

1 

i 

l 

5,DQPEN 

, CROSS 

,TIME 

,<<MA1  ,«X 

,«XS1 

,<XMIN 

,NDATF 

TRANP 

43 

5 

s 

\ 
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EPST=£PSILO*T0pP 
£RS7=.i 

xo=xp 
yo~yp 

ZO=ZP 
Tf>=TP 
SIGXO=0. 

STGYO=0. 

RO  =  0. 

NnftTC-NOftTr- 
Norci=o 
NDTOl=0 
KRHC1=0 
<BH01=0 
1000  CONTINUE 
<ay--i 

TP  *KRIP.£0.1>  GO  TC  50 

CALCULATE  FALL  VELOCITY  F CR  AOVECTION  ON  BASIS  OF  LOCAL 
C  FALL  RATE  AND  WIND  FIELD  AND  STORE  IN  WGAR 
CALL  CALI0<Z9H,KBKX,ZO,KAYfKBH:i 
CALL  CALI0CTIMUP,LTIMX,TO,1»LTIH) 
WBAR=WFZCICeHZ-i,N0AT0,LTIMl-CAYS««8HZ-i» 

IF  ( WBAP)  112,111,110 

CONFUTE  THE  VERTICAL  PSEUDO- VELOCITY  WBAR  ANO  STORE  ITS  SIGN  IN  RAY 
C  FOR  THE  OIFFUSIVE  SETTLING  CASE 
50  WOAR=(ZOEp-ZP) /(TDEF-TPI 
IF (WPARI  112,111,110 

110  XAY=<AY*1 

111  KAV=KAY*1 

112  CONTINUE 

CALIBRATION  OF  ZOEP  AGAINST  ZPH  YIELDS  TERMINAL  Z9H- PLANE 
IF  I <R IP.EQ.l)  CALL  CALTB<ZBH,*RHX ,Z0EP,<AY,<8H» 

CONSIOEK  KAY=0  CASE  INDEPENDENTLY 
200  IFfKAY.NE.O*  GO  TO  206 

IF  (RRIP.EQ.i)  GO  TC  205 

C  IN  THE  ADVECTIVE  TRANSPORT  CASE  WHENEVER  THE  ACTUAL  FALL  RATE 

C  IS  ZERO  THEN  SET  THE  DEPOSITION  TIME  INCREMENT  EQUAL  TO  THE 

C  TIME  LEFT  BEFORE  THE  WIND  FIELD  IS  UPDATED 
TSEG=TIMUP(LTIM*-l)-TO 
KBHC-K8HZ+1 
KBHO=KBHZ 
GO  TO  300 
205  TSEG-TOEP-TO 
MODE-1 
K3HC=K6H 
<BH0=28HC-1 
GO  TO  30C 

CALIBRATION  OF  ZO  AGAINST  ZBH  YIELDS  CURRENT  Z8H-PL ANf 
206  CALL  fyALIB  IZBh»  <BHX  ,iO,-)<AY,  KFRYI 
CALL  CAL'.BfZDH,KBHX,ZO,»<AY,<SHC: 

TF  (  i  XRIPvEQ.  0  > • OR  •  {KBHO.NE.  XTRYI  )  GC  .'0  2 5  3 
CONSIDER  EXCURSION  TO  TERMINAL  ZBH-PLANE 
iCBHC=<9u 
ZEST=ZEH»Ka>€) 

SOO£=fl 

1F£KAVMK8H0-'KBHC)  f IT  221,213,218 
CQN3IQE*  EXCURSION  BETWEEN  ADJACENT  ZDH-PlANES 
213  F6HC”XRHC>KAY 


TRANP  59 
TRANP  60 
TRANP  61 
TRANP  6? 
TRANP  6? 
TRANP  64 
TRANP  65 
TPANP  66 
TRANP  6? 
TRANP  6P 
TRANP  £9 
TRANP  70 
TRANP  71 
TRANP  7? 
TRAMP  72 
TRANP  74 
TPANP  75 
T°ANP  7f 
TRANP  77 
TRANP  73 
T "* ANP  79 
TR*'NP  30 
TRANP  81 
TPANP  82 
TRANP  83 
TRANP  84 
TRANP  85 
TRANP  36 
TRANP  97 
TRANP  88 
TRANP  89 
TRANP  SO 
TRANP  91 
TRANP  92 
TPANP  93 
TRANP  94 
TRANP  95 
TRANP  96 
TRANP  97 
TRANP  96 
7R»iMP  99 
TRANP1Q0 
TRANPiOt 
TRANP102 
TRANR103 
TRANP104 
TRANP105 
TRANP166 
TRANP107 
TRANP108 
TRAr:Pit)9 
TRANP110 
TRANPlil 
TRANPU2 
TRANP113 
TRANP114 
T° ANP115 
TRANPH6 
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UBAR, <PHF,NDATF,ITIMF 
VBAR,  KBMF,ND.ATF»LTIMF 


ZEST  =  Z:'H(<PHC) 

MOOE=l 

IF<  UBHC.,LE.<3HX>.  AKO.  {<AY*ZEST..LE.<AY*7DEP>)  GO  TO  221 
CONSIDER  EXCURSION  TO  TERMINAL  ZOEP-PLANE 
218  KBHC=<f?MO*<AY 
ZEST=ZOEP 
MOOC=1 

221  TSEG=<ZEST-ZO*/WBAR 
CHECK  IF  UPDATE  TIME  80UNCAPY  MILL  BE  CROSSED 
300  TC=TO*YSEG 

CALL  CALI«(TIMUP»LTINX,TO,l,LTIM» 

IF(  (LTIM.IT.LTIMX)  .  ANO.  (7IMUP  (LTIM*  i> .  LE.TCI  I  T  SEG=TlHUPC  LTIM+1 )  - 
COMPUTE  AVERAGE  HORIZONTAL  VELOCITIES  UBAR  AND  VBAR 
KDHA=KRHO 
<BHB=<RHC 

IF(KAY.LT.O)  GO  TO  40E 

<8HA=<PHC 

KBH8=K?H0 

405  CALL  GETOA  <  USUM,ZRF,K8HA, <BHB,NDATO,LTIM,  UBAR, <PHF, NDATF,ITIMF 
CALL  GETDAt  VSUH,Z«H,K9HA, KBHR , NDA TO,L TIM,  VBAR, KBMF,ND.ATF»LTIMF 
41i7  IF (NOEPjME.  0)  GO  TO  412 

COMPUTE  AVERAGE  HORIZONTAL  DISPERSION  AND  WIND  ORIENTATION  ANGLE 

CALL  C-ETDA  (OXSUM,  ZB*-,KRHA, KDHB,NDA TO, L f IM,DXBAR, XRHK, NDATF,LTIMr 
CALL  GETDA  (DYSUM,  Z8H?<DHA,  .<BHB,NGATO, lTIH ,OYE5AR,  <RKF, NDATF.LTIMF 
CALL  GETDA  (  RSUM,Z0F,KBHAsK8HR,NOATO;LTlM,  R8AR,  KRHF# NQATK,L  UMF 
RC=RO+RSAR 

SIGXC=SISXO*DXBAfc*TSEG 

SIGYC=SIGYO*OYBAR*TSEG 

COMPUTE  CURRENT  POSITION  ANO  TIME  (XC,YC ,ZC, TCi 
412  TC=TOfTSEG 

ZC=ZO+WBA?*TSEG 
XC  =  X(HU8AR*TSES 
yr  rVO*VBAR*TSEG 

CALL  NESTtNET#NETSU,XC,YC,NDATC,XL,XR,YL,YU.ICF, JCF,NCF| 

IF{MC(4) .EQ.l)  WPITE (ISOUT ,  31  TC,ZC,XC.  ,YCfNOATC 
COMPARE  CURRENT  MESH  INDEX  NDATC  WITH  PREVIOUS  MESH  INOEX  NDATO 
IF(NDATC.EO.NOATO)  C-C  TO  700 
COMPUTE  INTERPOLATED  POINT 
XT=XC 
YT  =  YC 
ZT  =  ZC 

IF  1SCDE.EQ.0)  GO  TC  21? 

CALL  aCUNfNET, wETSU,XT,YT,XO,YO,XC,YC, JCF,JCF,NCF1 
ZC=S0*7f  UYT-XC>**2  +  tYT-YC5**2»/UXT-XO)**2«-(YT-YO)*»2j* 
ZC=ZT+ZCMZO-ZTI 

IF<ApS(W8ARl.LE.l. 0E-3C)  GO  TO  510 


IF<ApS(WRARl.LE.l. TE-3C)  GO 
TSEG= (ZC-ZOJ/WBAR 
GO  TO  518 

IF' A6S(U8AR).LE.l. 0E-30)  GO 

TSEG=CXC-X0»7U8AR 

GO  TO  El  5 

TF<ABS<VBAfn.LE.U  PE-30?  GO 

TSEG~{YC-YO%/VBAR 

GO  TO  5i* 

CALL  ER?nRCPROGPN,5lE,ISOUn 
RETURN 

IF(NDEP.NE.O)  GO  TO  521 
RC=PO'R3AR 


GO  TO  5?  3 


TO  516 
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T»ANP119 
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)  TRAMP136 
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SIGXC=SIGXO*DXBAR*  T  SEG 
SIGYC  =  SIGYO«-DYBAP*TSEG 
921  TC-YO+TSEG 

CALL  NEST C  NFT  , NETSU  ,XC, YC»  NOATC  »XL  ,  XR,  YL»  YU*IC-F,  JCF»NCFI 
CHEC<  IE  PARCEL  CENTER  POSITION  IS  OSCILLATING 

TF ( (KQWOl.NE.KBHOJ .CR. (KBHC1.NE.K9HC1 .OR, (NDTCi. NE.NDATC1 .OP. 
1 (NDT01 »NE. NQATOl I  GO  TO  626 
IECMCt4l.EQ.il  WRITE  C ISOUT ,41  XP, YP , ZP , T» »XC. f C, ZC, TC 
CALL  CNTRt NET, NETSU  ,N0A70,XG,VG,ICE,JCF,NCF1 
XQ  =  XG 
Y  a  s  Y  6 

CALL  NEST  (NET,  NETSU  ,XQ,  YQ,  NGATQ,XLO ,  XRO,Yt.O,YUO,  ICF  ,JCF,NUFJ 
CLEAR  SYOREO  MESH  AND  STRATUM  INDICES 
NnTCl=0 
NO  TO 1- 0 
<OHCl=7 
KOHOi-O 

channel  wafer  center  position  along  aopppoiate  3Ell  bounce* 

SPE=2.*rPS 

IFMABS(XLO-XR).GT.SPE).»ANO.  (ABSCXPO-XLI.GT.S°f>  J  GO  TO  616 

UOAR-'O, 

CALL  GETOA  t  VSUM,ZOH,KRFA,  KBHP,fJOATO,LTTM,VBdPC,  KB.'ir  ,\'9ATF,Lr 

IF  ( A^S  (V8ARC)  .  LE .  A3?  \  V.'ARt  I  GO  TO  407 

VKaR=VRARC 

NOATO=NOATC 

GO  TC  407 

616  IEmaBS(YLO-Y:<).GT.SPE).ANO.  <  ABSO  UQ-  rL >  . GT.SPE)  I 
1CALL  ERROR (PROGPM, 626 , ISOUT1 
VPAR-sQ. 

CALL  GETOA  (  US!J*,Z9f,  K9HH,N0A7G,LTTH»UeARC,  RBHF,  N!?ATF,l  \ 

IFCABSCUBALC)  .LE.ABE(UEA»n  GO  TO  407 

UBrtP=UBARC 

fjQATC=NOATC 

GO  TO  407 

COMMIT  PREVIOUS  ANO  CUR^tHT  MF3H  ANC  STRATUM  INDICES  TO  STORAGE 
626  NOTC1=NOATO 

notoi=noltc 

KBr:C  1  “XPHC 
KBHCl=<PHO 

CONVE°T  XO,YO,ZO, TO,SIGXO,S1GYO,  AND  NO A  TC  TO  CU^nT  V-LUES 
700  ZO=ZC 

ro=xc 

YO=YO 

TO=TC 

noato=ndatc 

IF  CMC (41 tEO.13  WRI TE ( i SOUT ,21  TO, Z C, XO,YO,NOATO 

IF(NOEPnNE.O)  GO  TC  70S 

SIGXG-SIGXC 

SIGY0=SIGYC 

RO=RC 

CHECK  IE  CURRENT  POSITION  IS  OUTSIOE  4TM0SPHERE 
708  IF (NOATU.LE.G?  GO  TO  710 
IE  (KRIP.EQ.t)  50  T C  7C9 

C  If  DEPOSITION  PLANE  XS  REVOKED  OR  TRANSPORT  Tlftf  LIMIT  IS  €X» 

C  EXIT  FROM  TRANP,  OTHERWISE  RETURN  TO  » OP 

IE  J(  (70-7DEPl  .LE.EPS'21  .OR.  *  ("IMSX- fG* . LE.EPSTI  i  GO  TO 

GO  TO  10AG 
709  CONTINUE 


-RANP17S 
T°ANPl76 
TRANP177 
TEANPl'B 
TPANP179 
TPANP130 
TOANP131 
TRANP182 
TRAMP183 
TRANP184 
Y°ANP185 
TRANP1B6 
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TPANP202 
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IMf)  TRAMP204 
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TOANP212 
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TRAMPS J1 
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CHECK  IF  TOTAL  FLIGHT  TIME  HAS  REEN  EXHAUSTED 
IF (TO*EPST«TOEpt  300,720,720 

CARRY  PAOCEL  CENTER  RACK  INTO  AT40S.  IF  NDEP  IS  NOT  EQUAL  TO  ZERO 
710  IF (NOEP.EQ.  01  GO  TO  720  '  >' 

XO=XO-2.*EPS. 

YO=Y0-2.*tPS 

CALL  NEST (NET, NETSU ,  XO , YO, NOATO, XL , YP, YL, YU,  ICF, JCF,NCF) 
IFCNOATO.GT.O)  GO  TC  720 
714  XO=XO«-A.  «EPS 
Y0=Y0+4.*EPS 

CALL  NEST?  NEr , NETSU, XO,YO«  NO A TO, XL ,TR»YL»YU,ICF, JCF,NCF» 

IF { NnftTC.LE.OJ  CALL  £PROR(PROGPM,720,ISOUT) 

COMPUTE  HOP I  Z«  DISPERSION  IF  NOCP  IS  NOT  EQUAL  TO  ZERO 

720  IF{NOEP.NF,0>  RETURN 
P2=PWAP**2 

IFtMCdOi.EQ.l)  GO  TO  721 
Sinxc=2* ^OGKN*?IGXC 
STGYC=2.*CROSS*SIpYC 
GO  TO  722 

721  rprp=rc-rp 

OSPPTX=SIGXO/T  RIP 

rONEX  =  4.M«R2'OSPRrX)*«U,/?.n 

Ir  (TRIP.LE.TONEXI  Slr*xC=0SPRTX*T0NEX*(TRIF**2)  73. 

IF (TRIP. GT. TONEX)  SIGXC=0SoPrx7(rRIP*^3l/3. 
STGXO=SIGXO*(DOHN**  (3./2.1 > 

OSPRTY=SIGYOATRIP 

TONEY=4. *(  CR270SPRTY)**(1.73.)) 

IF ( TRIP. LE. TONEY)  SIGYC=OSPRTY*TONjY*<TPfP**2»  '3. 
IF(TPIP.GT. TONEY)  S 1GYC-DSPRTY* (TRI°**3V 73. 
SIGYO=SIGYO«(CROSS**<3./2.>) 

722  SIGX0=S0RT(R2«-SIGX0I 
SIGYO=SORT(R2^SIGYOI 
RETURN 

ENfl 
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